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FINITE  ELEMENT  MODELING  OF  STATIC  INDENTATION  DAMAGE  IN 
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BY 
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Chairman  :  Bhavani  V.  Sankar 
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The  objective  of  this  research  is  to  study  progressive 

failure  of  circular  laminated  brittle  matrix  composite  plates 

due  to  static  indentation  loading  using  finite  elements.  The 

plates  are  assumed  to  be  simply  supported  along  the  edge,  and 

are  indented  by  a  rigid  indenter  with  a  hemispherical  nose.  A 

quasi-isotropic  homogenization  process  is  used  to  make  the 

problem    axisymmetric.     Four    noded,     isoparametric, 

quadrilateral,  axisymmetric  solid  elements  are  used  to  model 

the  plate.  A  new  algorithm  is  used  to  analyze  the  progressive 

contact  between  the  indenter  and  the  plate.  Hashin's  failure 

criteria  are  modified  so  that  they  can  be  applied  to 

axisymmetric,  transversely  isotropic  laminates.  Stiffness 

matrices  of  failed  elements  are  modified  suitably.  New  slip 

surface  is  introduced  when  delamination  occurs.  Very  stiff 
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springs  are  attached  to  nodes  on  delamination  surfaces  to 
prevent  interpenetration  of  the  nodes.  The  results  from 
incremental  damage  analysis  agree  well  with  experimental 
results.  This  program  can  be  applied  to  analyze  impact 
problems  involving  large  masses  impacting  at  very  low 
velocities. 
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CHAPTER  1 
INTRODUCTION 

1.1   Introduction 

Currently  composite  materials  are  used  in  primary- 
structures  of  aircraft,  particularly  aircraft  wings  and 
frames,  because  significant  weight  savings  can  be  achieved  due 
to  their  high  strength  to  weight  ratio.  Also,  composite 
materials  can  easily  be  tailored  to  obtain  desired  properties 
in  different  directions.  For  example,  more  0°  layers  can  be 
used  for  longitudinal  loads,  90°  layers  for  transverse  loads 
and  ±  45°  layers  for  shear  loads.  Also,  there  are  other 
advantages  such  as  resistance  to  damage  by  fatigue  loading, 
immunity  to  corrosion,  ease  of  processing,  and  low  coefficient 
of  thermal  expansion.  Structural  forms  that  are  otherwise 
inconvenient  or  impossible  to  manufacture,  e.g.,  more  complex 
airfoils  or  helicopter  blades,  can  be  easily  fabricated. 
Hence,  the  utilization  of  composite  materials  in  aerospace  and 
transportation  industries  is  continuously  increasing  [1,2]. 

However,  it  is  well  known  that  brittle  matrix  composites 
such  as  graphite  epoxy  are  very  susceptible  to  low  velocity 
impact,  which  may  be  caused  by  dropping  tools,  runway  debris, 
hailstones,  bird  strikes  and  ground  service  vehicles  hitting 
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aircraft  structures.  High  and  medium  levels  of  impact  energy 
cause  surface  damage  that  is  relatively  easily  detected,  and 
therefore,  repairs  may  be  undertaken.  However,  low  velocity 
impacts  which  are  at  quite  low  energy  levels  can  seriously 
damage  laminated  composite  structures  internally  so  that 
frequently  it  is  very  hard  to  detect  damage  through  visual 
inspections.  The  damage  can  be  detected  only  by  using 
nondestructive  inspection  (NDI)  techniques.  For  example,  low 
velocity  impact  to  an  aluminum  plate  may  cause  denting  damage 
and  the  structure  does  not  lose  very  much  strength.  However, 
inside  of  an  impacted  composite  plate  there  may  be  fiber 
breakage,  matrix  cracks,  and  delaminations,  while  the  impacted 
surface  may  not  show  any  damage.  Commonly  it  is  termed  as 
"Barely  Visible  Impact  Damage"  or  BVID  [2].  After  impact, 
composite  materials  lose  considerable  strength  because  of 
internal  damage  and  delaminations.  [2-7]  >> 

It  becomes  very  important  to  understand  the  response  of 
composites  to  the  impact  of  rigid  objects  and  to  predict 
damage  in  order  to  make  composite  materials  more  reliable.  The    '  . 
study  of  the  effect  of  BVID  on  static  and  fatigue  strengths  of 
graphite/ epoxy  panels  shows  that  impact  degradation  is  of  more 
concern  in  compression  loading  than  in  tension  because  of   "* 
delaminations  which  occur  during  impact.  Importance  of  the      jj 
impact  study  is  emphasized  by  the  fact  that  upper  surfaces  of 
the  aircraft  are  more  likely  to  receive  impacts  from  tool 
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drops  than  lower  surfaces,  and  upper  surfaces  operate  in  a 
compression-dominated  load  regime  [2]. 

Because  of  the  extremely  complicated  nature  of  impact  on 
the  composite  materials,  empirical  methods  have  been  widely 
employed.  Generally,  empirical  methods  have  several  distinct 
limitations.  First,  it  is  very  expensive  and  also  time 
consuming  to  conduct  testing  and  to  process  the  data.  Second, 
it  is  practically  difficult  to  test  all  the  combinations  of 
material  properties,  boundary  conditions  of  the  composite 
structures,  and  impactor  velocities  and  weight  variations. 
This  problem  can  be  critical,  if  the  composite  material  is  in 
a  developing  stage,  and  only  limited  quantities  are  available. 
Third,  it  is  very  difficult  to  comprehend  all  the  stress  and 
strain  distributions,  and  damage  progression  inside  composites 
by  empirical  methods.  And  also  there  are  difficulties  in 
matching  results  from  specimens  to  the  real  structures.  [6,7] 
Therefore,  numerous  analytical  methods  have  been 
developed  to  overcome  the  shortcomings  of  empirical  methods. 
However,  so  far,  because  of  the  extremely  complicated  nature 
of  impact  problem,  there  are  no  complete  three-dimensional 
models  available  which  can  simulate  the  entire  impact  process 
and  predict  damage  progression.  There  is  an  urgent  need  to 
develop  a  three  dimensional  progressive  finite  element  program 
which  calculates  contact  forces,  contact  displacements, 
stresses  and  strains  inside  the  composite  plate.  During  the 
contact  process,  the  finite  element  program  must  be  able  to 


4 
check  all  the  failure  criteria  for  fiber  breakage,  matrix 
cracking  and  delamination.  If  failure  occurs,  then  the  program 
could  modify  the  stiffness  matrix  to  account  for  failures  . 

Fortunately,  low  velocity  impact  due  to  large  masses  can 
be  treated  as  a  static  indentation  problem,  because  the  impact 
duration  is  much  longer  than  the  time  required  by  the 
propagating  waves  to  travel  from  the  impact  site  to  the 
supports  or  free  edges,  and  the  test  results  show  good 
agreement  with  this  assumption.  Strictly  speaking,  the 
terminology  of  low  velocity  impact  is  not  clearly  defined. 
Rather  low  impact  velocity  is  generally  defined  as  the  impact 
velocity  which  is  low  enough  to  justify  a  static  analysis  of 
the  response  of  the  structure,  especially  in  the  vicinity  of 
impact  [3,4,8,9].  This  velocity  depends  upon  the  mass,  shape 
and  the  material  property  of  the  impactor  and  also  upon  the 
impacted  structure  boundary  conditions  and  its  natural 
frequencies.  Generally  it  ranges  1-3  meter/second  with 
impactor  mass  less  than  10  kg  [10]. 

The  advantage  of  static  analysis  is  that  the  time 
variable  is  eliminated.  Hence,  we  can  closely  look  into 
progressive  damage  mechanisms  with  a  small  amount  of 
computational  resources.  If  time  effects  are  considered,  it 
will  be  a  formidable  task  to  obtain  all  the  meaningful  stress 
and  strain  distributions  inside  the  plate  and  to  map  the 
damage  pattern  and  its  progression  because  of  computational 
time  and  storage  limitations  of  the  computer  system.  If 


5 
coarser  elements  are  employed,  it  will  be  very  hard  to  obtain 
accurate  data  because  meaningful  data  must  be  obtained  in  the 
small  confined  area  under  the  indenter. 

The  experimentally  determined  contact  laws  account  for 
local  damage  in  an  empirical  manner,  and  hence  cannot  be  used 
to  study  the  local  damage.  A  fully  three  dimensional 
formulation  is  required  to  study  the  progressive  damage. 

In  the  present  study  a  finite  element  program  performs  an 
incremental  analysis  of  the  contact  problem,  and  assesses  the 
damage  due  to  indentation. 

1.2  Experimental  Background 

In  this  section  some  results  from  an  experimental  study 
[11]  are  presented  to  bring  out  the  importance  and  necessity 
of  the  finite  element  modeling  of  the  damage  process  due  to 
static  indentation. 

Test  specimens  were  made  of  Hercules  carbon  prepreg  tape 
AS4/3501-6,  which  is  an  amine-cured  epoxy  resin  reinforced 
with  unidirectional  carbon  fibers.  All  specimens  were  cut  from 
a  30.48  cm  (12  inch)  by  30.48  cm,  [0°,  +45°,  -45°,  90°]35 
composite  plate  lay-up.  In  order  to  obtain  the  specimens  with 
the  quality  suggested  by  the  manufacturer,  the  curing 
instructions  was  carefully  followed  as  in  Ref.  [12].  An 
autoclave  (Baron  Model  BAC-24)  was  used  to  fabricate  all  the 
specimens  of  composite  laminates  used  in  this  study.  Detailed 
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procedure  of  fabrication  can  be  found  in  Ref.  [11].  Material 
property  of  the  0°  lamina  can  be  found  in  Table  2.1. 

Figure  1.1  shows  the  pendulum  impact  test  facility  in 
which  all  impact  tests  were  conducted.  The  pendulum  includes 
the  support  strings  connecting  the  tup  to  the  laboratory 
ceiling,  and  the  tup,  which  is  a  combination  of  weights,  a 
load  cell,  and  an  indenter  with  one  inch  diameter 
hemispherical  head.  The  total  mass  of  the  tup  is  13.84  kg. 
Static  indentation  tests  were  done  in  a  standard  MTS  material 
testing  machine. 

The  composite  laminate  was  placed  on  a  circular  steel 
ring  of  5.08  cm  (2  inch)  diameter  to  simulate  a  circular  plate 
with  simply  supported  boundary  condition.  The  test 
configuration  is  shown  in  Figure  1.2. 

The  test  results  of  contact  force  vs.  center  deflection 
of  the  circular  composite  plate  for  the  static  indentation  and 
impact  tests  are  shown  in  Figure  1.3.  This  graph  suggests  that 
static  indentation  tests  can  simulate  low  velocity  impact 
tests.  Also  it  suggests  that  indenting  process  is  highly 
nonlinear  due  to  internal  damages.  A  sudden  drop  of  the 
contact  force  is  a  unique  feature  of  indentation  of  composite 
plates.  Figure  1.4  shows  the  picture  of  the  C-scan  of  the 
impacted  plate.  The  circular  damage  pattern  justifies  the 
axisymmetric  formulation  used  in  the  present  study.  Figure  1.5 
shows  the  picture  of  the  photo-micrograph  of  the  polished 
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cross  section  of  the  indented  plate.  This  shows  matrix  cracks 
due  to  shear  failure  and  delaminations. 

However,  these  data  are  not  sufficient  to  understand  the 
entire  impact  or  indentation  process.  They  do  not  explain  how 
damage  initiates  and  propagates.  Also,  it  is  very  hard  to 
quantify  this  damage. 

1.3   Objective  and  Scope 

The  major  objective  of  this  research  is  to  develop  a 
finite  element  program  which  models  the  static  response  of 
laminated  composite  plates  indented  by  a  rigid  spherical 
indenter  and  to  predict  progression  of  damages  inside  the 
plate.  The  results  are  applicable  to  low  velocity  impact  of 
the  composite  plates.  '    r v 

A  circular  plate  and  a  rigid  hemispherical  indenter  were 
chosen  in  modeling  because  they  correspond  to  many  practical 
impact  situations.  Quadrilateral,  isoparametric,  axisymmetric 
finite  elements  were  employed  in  this  research.  By  employing 
axisymmetric  elements,  the  three  dimensional  problem  can  be 
treated  as  two  dimensional  by  eliminating  dependency  on  6. 

In  Chapter  2,  the  finite  element  fomnulation  and  the 
incremental  displacement  method  are  introduced.  The  input  to 
the  finite  element  program  is  the  displacements  necessary  to 
close  the  gap  between  the  indenter  and  candidate  contact  nodes 
of  the  plate.  The  laminated  plate  consists  of  many  plies  that 
are   orthotropic   and   have   different   angles   of   fiber 
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orientation.  As  axisynmetric  finite  elements  were  employed  in 
this  research,  there  must  be  special  considerations  of  the 
material  properties  and  failure  criteria  to  account  for  this 
discrepancy.   The  in-plane  strains,   e„,   e^  and  Yw.  ^^s 

XX      yy  xy 

continuous  at  the  interface  between  two  adjacent  plies  whereas 
the  stresses,  o^, ,  o^  and  x^^  are  not.  On  the  other  hand,  o,,, 
T^^  and  T  are  continuous  through  thickness  and  are  better 
represented  in  a  homogeneous  plate  model  than  e^^,  y,^^  and  y  , 
which  can  be  discontinuous  through  the  thickness.  A  set  of 
transversely  isotropic  elastic  constants  is  derived  for  a 
group  of  plies  that  can  be  considered  as  a  quasi-isotropic 
laminate. 

In  Chapter  3 ,  the  accuracy  and  efficiency  of  the  present 
method  were  evaluated  by  comparing  the  results  with  some  early 
work  [10].  The  contact  force  vs.  center  deflection  curve  of  a 
AS4/3501-6  graphite/ epoxy  plate  was  compared  with  experimental 
results  in  the  elastic  range.  Then,  some  widely  used  contact 
laws  were  compared  with  the  present  finite  element  results. 
The  effect  of  friction  was  studied  by  assuming  stick  friction 
conditions  between  the  indenter  and  the  top  surface  of  the 
plate. 

In  Chapter  4,  the  failure  criteria  are  introduced.  Among 
many  interactive  and  independent  failure  modes,  three- 
dimensional  failure  criteria  introduced  by  Hashin  [13,14]  were 
employed  in  the  present  studies.  These  criteria  have  the  same 
quadratic   form  as   Tsai-Hill's   criterion,   but  they  can 
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distinguish  different  failure  modes  such  as  fiber  breakage, 
matrix  cracking  and  delamination.  Since  we  use  an  axisymmetric 
formulation,  the  damage  pattern  is  also  going  to  be 
axisymmetric,  which  does  not  strictly  represent  the  damage 
pattern  in  a  composite  laminate.  However,  experimental  impact 
studies  [15]  have  indicated  that  the  actual  damage  pattern  is 
like  a  spiralling  stair  case,  and  the  projection  of  these 
delaminations  can  be  considered  as  circular  and  hence 
axisymmetric.  This  was  also  confirmed  by  studies  in  the 
laboratory  [11].  The  extent  of  delamination,  fiber  breakage 
and  matrix  failure  were  mapped  during  the  contact  process. 
After  failures  are  detected,  the  total  stiffness  matrix  of  the 
plate  must  be  modified  properly  depending  upon  the  failure 
modes.  Matrix  and  fiber  failure  can  be  implemented  by  reducing 
or  eliminating  element  stiffness  matrices.  Implementing 
delamination  into  the  finite  element  program  requires  creation 
of  new  free  surfaces,  consequently  new  nodes.  Rigid  springs 
are  inserted  across  the  newly  created  free  surfaces  to  prevent 
interpenetration  of  the  nodes.  Delamination  can  be  assumed  as 
Mode  II  (Shear  Mode)  crack  propagation.  Incremental  damage 
analysis  results  are  compared  with  experimental  results. 

A  summary  of  conclusions  and  some  suggestions  for  future 
research  are  listed  in  Chapter  5. 
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Figure  1.1   Impact  test  Scheme 
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Figure  1.2   Test  configuration 
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Figure  1.3  Contact  Force  vs.  center  deflection  Curve 
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Figure  1.4   C-Scan  picture  of  impacted  composite  plate 
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Figure  1.5  A  photo-micrograph  of  the  impacted  plate 


CHAPTER  2 
FINITE  ELEMENT  PROGRAMMING  OF  THE  CONTACT  PROCESS 

2 . 1   Introduction 

The  finite  element  method  has  been  established  as  a 
powerful  and  versatile  tool  of  analysis  for  contact  problems 
involving  isotropic  bodies.  In  the  present  study,  a  method  for 
solving  the  problem  of  indentation  of  composite  laminates  is 
developed. 

The  contact  problem  of  quadratic  surfaces  of  elastic 
bodies  was  first  solved  analytically  by  Hertz  in  1881.  Since 
then,  numerous  contact  problems  have  been  solved  by 
formulating  the  problem  as  an  integral  equation.  The  solution 
is  either  obtained  by  solving  the  equation  numerically  or,  for 
simple  configurations,  in  closed  form.  The  Hertzian  type 
problem  has  three  general  assumptions.  First,  the  contacting 
solids  are  homogeneous,  isotropic  and  linear  elastic.  Second, 
the  contact  areas  are  essentially  flat  and  relatively  small 
compared  to  the  various  radii  of  curvature  of  the  undeformed 
bodies  in  the  vicinity  of  the  contact  interfaces.  Third,  the 
contacting  solids  are  perfectly  smooth  and  only  normal 
pressures  are  considered  [6,16]. 
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Hertzian  type  solutions  are  not  suitable  for  indentation 
of  composite  laminates  which  are  inhomogeneous,  anisotropic 
and  also  involve  significant  damages  at  low  load  levels. 
Further,  most  laminated  composites  in  use  can  not  be 
adequately  represented  as  a  half-space  as  in  Hertzian  contact 
problem.  Because  of  these  complexities,  three  dimensional 
analysis,  especially  three  dimensional  finite  element  methods, 
are  required  in  the  present  study. 

In  the  last  two  decades,  many  researchers  thought  that 
contact  problems  could  be  considered  as  a  special  case  of 
constrained  minimization  of  either  total  or  complementary 
potential  energy.  The  minimization  is  formulated  as  a 
mathematical  programming  problem,  and  the  solutions  are 
obtained  by  either  using  incremental  linear  programming, 
quadratic  programming  or  conjugate  (modified)  gradient 
methods.  The  minimization  problem  can  also  be  formulated  as  a 
variational  inequality  using  penalty  methods  or  Lagrange 
multipliers  [17-21].  However,  these  methods  have  certain 
disadvantages  such  as  creation  of  flexibility  matrix,  the 
awkwardness  involved  in  introducing  additional  unknowns,  such 
as  slack  variables  and  Lagrange  parameters.  Nothing  can  be 
better  than  direct  formulation  of  finite  element  method  as  a 
simple  design  tool,  which  calculates  the  detailed  stress  field 
in  various  types  of  laminates  for  various  boundary  conditions 
accounting  for  the  damages  inside  the  plate.  The  goal  is  to 
use  simple  and  conventional  elements  so  that  the  method  can  be 
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implemented  in  commercial  finite  element  programs.  The  novel 
incremental  contact  algorithm  developed  in  this  study  is 
explained  in  Section  2.4. 

2.2  Axisymmetric  Formulation  - 

The  benefit  of  employing  an  axisymmetric  element  is  that 
the  three  dimensional  problem  can  be  treated  as  two 
dimensional.  Especially  in  the  case  of  rigid  ball  indentation, 
where  the  contact  surface  is  circular,  axisymmetric 
formulation  is  a  very  efficient  method  to  treat  the  indenting 
process  in  a  computer  program,  because  upper  surface  nodes  of 
the  plate  come  into  contact  successively  and  the  contact  area 
is  readily  determined.  If  we  use  three  dimensional  solid 
elements,  it  will  be  a  formidable  task  to  determine  the 
contact  nodes  for  the  case  of  a  spherical  indenter,  although 
the  employment  of  three  dimensional  solid  elements  has  a 
distinct  advantage  in  that  we  can  use  orthotropic  material 
properties  and  also  the  failure  criteria  for  orthotropic 
laminates  directly  without  any  modifications.  Complexity  of 
the  contact  algorithm  demands  employment  of  axisymmetric 
elements,  and  consequently  there  must  be  some  considerations 
to  convert  the  material  properties  of  a  laminated  plate,  which 
are  basically  orthotropic,  into  equivalent  axisymmetric 
properties.  The  circular  damage  distribution  pattern  of  the  C- 
scan  pictures  of  the  damaged  plates  (Figure  1.4)  justifies  the 
employment  of  axisymmetric  elements  to  solve  the  problem. 
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Stress-strain  relationships  of  the  axisymmetric  element  are 
discussed  extensively  in  Section  2.3. 

In  the  present  study  u(r,z)  and  w(r,z)  represent  the 
radial  and  transverse  displacements,  respectively.  The 
axisymmetric  stress  and  strain  components  are 


(2.1) 


{   a   }    =   {   a^^     a^^    t^^  o^q     } '' 
(e}  ={  e,,  e,,    Yrz  ^ee  }^  (2.2) 


where. 


^rr 

du 
dr 

^zz 

dw 
dz 

^TZ     = 

du   , 

dz 

dw 
dr 

^ee 

u 
r 

(2.3) 


As  graphite/epoxy  composite  is  very  brittle  and  can  be 
considered  as  elastic  until  failure,  material  non-linearities 
are  ignored  [2].  As  will  be  seen  later,  the  deflections 
considered  are  smaller  than  the  plate  thickness,  and  hence 
geometric  non-linearities  are  also  not  considered. 
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2.3  Constitutive  Relations  for  an  Axisymmetric  Laminated  Plate 

A  composite  laminate  can  be  considered  as  having  several 
sublaminates.  If  each  sublaminate  is  quasi-isotropic,  then 
each  sublaminate  can  be  idealized  as  a  homogeneous 
transversely  isotropic  material.  Then  the  composite  laminate 
can  be  considered  as  a  laminate  made  up  of  several 
transversely  isotropic  layers.  Hence  the  axisymmetric 
formulation  can  be  employed  if  the  external  loading  and 
boundary  conditions  are  also  axisymmetric. 

The  stress-strain  relationship  in  the  Cartesian 
coordinates  are 

{e  }  =  [  5]  {  a  }  (2.4) 

where  (  o  }  and  {  e  }  are, 

{  e  }  =  {  e^  63^  Y^  e^z  Y^x  Yzy  }''  (2.5) 

{  O  }  =  {  ^xx  «^yy  "^xy  <^ zz    ^ zx    "^ zy  V  ^^-S) 

For  the  ply  which  has  fiber  orientation  angle  6  from  the 
X  coordinate,  lamina  constitutive  relations  can  be  obtained 
from 

{ejo  =[  7;]^[5][  r„  ]  {  a  }e  (2.7) 


where    [    S    ]    is 
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and    [    T     1  ,     [    T     ]    are   given  by, 
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For  the  fiber  orientation  angle  9, 


{  e  U  =  [  S  ]q  {  o  )q 


[S],=[T,],'[S][T^]^ 


Averaging  [  S  ]g  for  the  number  of  plies  from  this  form 
is  not  suitable  because  in-plane  stresses  such  as  o^^,  a^  and 
X        are  not  continuous  through  thickness  of  the  element  while 

xy 

e  ,   G    and  y    ^^^  continuous .   And ,   also ,   out-of -plane 

XX     yy  xy 

stresses  o  ,  t  and  -c  are  continuous  through  thickness  while 
out-of-plane  strains  e^^,  y^x  ^^^  Y^y  ^^^  ^o"*^  continuous.  These 
are  well  illustrated  in  Figure  2.1  and  Figure  2.2,  which  Wu 
and  Springer  have  obtained  [23].  These  figures  show  the 
continuity  of  the  stresses  and  strains  through  the  thickness 
of  the  impacted  plate.  The  best  way  to  average  the  modulus  of 
the  element  can  be  done  after  denoting  [  S  ]q  as. 


[•5]e 


P 


(2.8) 


where  p,  q,  r  and  s  are  3x3  matrices.  Then  the  stress-strain 
relationship  can  be  written  as 
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The  stress-strain  relationship  can  be  rearranged  such  that  all 
the  components  of  stress  and  strain  that  are  continuous 
through  the  thickness  are  on  one  side  and  the  remaining 
components  are  on  the  other  side  as  shown  in  Figure  2.9. 
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(2.10) 


Let  us  denote  [  D  Jg^ 


[^]e 


P 


-1 


-p'^g 


rp  ^  s-rp~'^g 


(2.11) 


Then  an  average  [  D  ]  matrix  for  a  group  of  N  plies  can 
be  obtained  as, 
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•^   i=l 


P  Q 

R   S 


(2.12) 


Now  the  elastic  matrix  [    C    ]^^   can  be  recovered  from 
[  D  ]  .   The  [  C  ]   can  be  expressed  in  terms  of  P,  Q,  R  and 
S  as, 


C]e    = 


(2.13) 


The  stress-strain  relations  take  the  form 

For  widely  used  cross  ply  or  quasi-isotropic  laminates, 
such  as  [  0  ,  90  Is  or  [  0°  ,±  45°  ,  90°  Ig  laminates,  [  C  j^^ 
becomes  transversely  isotropic.  For  axisymmetric  problems, 
where 


^re  =  o^  =  e^Q  =  e^e  =  0 


(2.14) 
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[  C],  = 
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(2.15) 


Further,  the  present  homogenization  scheme  yielded  a 
symmetric  stiffness  matrix,  i.e.,  C.-^   =  C^.,    and  also        .  f 


(2.16) 


Sometimes,  five  engineering  elastic  constants,  such  as 


E,,  E^ 
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'   can  describe  the  material  properties  of  a 


transversely  isotropic  laminate  plate  [10].  For  this  case,  the 
stress-strain  relationship  can  be  obtained  as  follows. 
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Table  2.1  shows  the  material  properties  of  the  AS4/3501-6 
graphite/epoxy  lamina.  Also,  the  calculated  transversely 
isotropic  element  material  property  matrix  [  C  j^^  for  a  plate, 
which  consists  of  four  plies  in  [  0°  ,  ±  45°  ,  90°  ]  orientation 
is  given  in  Table  2.1.  This  configuration  is  chosen  because  of 
its  wide  application  in  aircraft  structures.  For  example,  the 
main  wing  box  for  a  straight  wing  aircraft,  0°  plies  carry  the 
spanwise  direct  stresses  induced  by  wing  bending  like  a 
cantilever  beam  bending,  ±  45°  plies  carry  shear  stresses 
caused  by  wing  twisting,  and  finally  90°  plies  carry  the 
chordwise  direct  stresses  due  to  fore  and  aft  bending  of  the 
aircraft  wing  [2]. 

2.4  Algorithm  for  the  Contact  Problem 
Unlike  conventional  finite  element  problems,  contact 
problems  are  mixed  boundary  value  problems  where  nonzero 
displacements  are  specified  at  parts  of  the  boundary  and  zero 
tractions  are  specified  in  the  rest.  Further,  the 
displacements  in  the  contact  surface  have  to  satisfy  certain 
conditions  specified  by  the  geometry  of  the  indenter.  Sankar 
[10]  used  a  contact  algorithm  in  conjunction  with  a  finite 
difference  method.  In  the  present  study,  the  algorithm  has 
been  modified  to  suit  the  finite  element  method.  Unlike  the 
finite  difference  method,  the  finite  element  method  has  the 
advantage  of  modeling  complex  geometries  and  multi-layer 
media. 
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For  the  frictionless  contact  problem,  the  vertical 
degrees  of  freedom  of  the  top  surface  nodes  of  the  plate 
become  contact  candidates.  In  Equation  (2.18),  the 
displacements  of  these  contact  candidate  degrees  of  freedom 
are  notated  as  d^^  and  other  displacements  as  d^.  The  d^  degrees 
of  freedom  can  be  condensed  to  obtain  a  smaller  stiffness 
matrix  as  shown  below. 


^u 

K^ 


0 


(2.18) 


{   d,}    =   -[  K^,]-^    {  K^^]     {  d„  } 


(2.19) 


After  condensation,  the  equilibrium  equations  take  the  form 


[  K^]    {  d„  }  =  {  F^  } 


(2.20) 


where 


[  i^.,  ]  =  [  ic„„  ]  -  [  ic„^  ]  [  i<r„  ]  -M  K^^  ] 


(2.21) 


Since  we  are  dealing  with  an  axisymmetric  problem,  the 
contact  area  is  always  circular.  Let  us  assume  that  the 
contact  has  progressed  to  Node  k-1  (see  Figure  2.4).  The 
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procedure  for  determining  the  incremental  load  required  to 
bring  the  k"^^  Node  into  contact  is  as  follows.  Equation  (2.20) 
can  be  partitioned  as  -■  W  ;  . 


W^ 


w. 


(2.22) 


where  w  are  the  z  displacements  of  nodes  within  the  contact 

c 

region  including  k-1  and  w^  are  the  remaining  degrees  of 
freedom  for  which  the  forces  are  zero.  Let  us  use  g,.  to  denote 
the  gap  between  the  indenter  and  nodes  outside  the  contact 
region,  i.e.,  i  ^  k.  The  first  step  is  to  determine  the 
incremental  displacement  of  the  indenter  to  close  the  gap  g,^. 
We  first  apply  a  unit  displacement  to  the  indenter,  i.e., 


Wj^  =  l,i  =  l,2,...,k-1 


(2.23) 


The  equation  (2.22)  is  solved  for  F^.  and  w^.  Now  the  new  gap 
q  • .  becomes 


g'.=g^-l+w^,     i  h  k 


(2.24) 


,  t 


Now  we  can  calculate  the  displacement  6  of  the  indenter 
to  make  the  gap  g',^  identically  equal  to  zero. 


'it-  !■  ■> 
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6  =     ^^ 


"  (2.25) 

9k 


(  1  -  V.  ) 


After  determining  6,  the  nodal  force  F^  can  be  computed.  By 
integrating  all  the  nodal  forces,  we  obtain  the  incremental 
contact  force,  AP.  The  procedure  is  repeated  for  the  next 
contact  node  k+1. 

When  frictional  effects  are  included,  radial  forces  at 
nodes  in  the  contact  region  can  not  be  considered  as  zero.  The 
procedure  has  to  be  slightly  modified,  as  is  explained  in 
Section  3.5. 

The  damage  alters  the  stiffness  matrix.  We  assume  that 
the  indentation  process  is  under  displacement  control,  and 
hence  the  w-displacements  of  the  contact  nodes  are  kept  the 
same  when  passing  from  one  iteration  to  the  next.  The  changes 
in  [  K  ]  cause  a  change  in  the  contact  load  (usually 
reduction)  for  the  same  indentation  displacement.  Detailed 
damage  analyses  are  presented  in  Chapter  4 . 

2.5   Axisymmetric  Finite  Element 
In  this  section,  the  derivation  of  the  stiffness  matrix 
of  a  quadrilateral,  isoparametric,  axisymmetric  element  is 
given  for  the  sake  of  completeness. 
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In  isoparametric  elements,  the  element  coordinates  and 
the  element  displacements  are  expressed  with  the  same  shape 
functions  associated  with  nodal  coordinates  and  displacements 
respectively  as  shown  below. 


i=l 
n 

1=1 
n 

i=i 
n 


(2.26) 


2  =  1 


where  r  and  z  are  element  coordinates,  ^  and  ti  are  the  natural 
coordinates  and  u  and  w  are  element  displacements  as  shown  in 
Figure  2.5.  Subscript  i  signifies  node  number.  Shape  function 
or  interpolation  function  N^  are  expressed  in  the  natural 
coordinate  system,  which  has  variable  ^  and  tj  that  each  vary 
from  -1  to  +1  as  shown  in  Figure  2.5.  The  characteristic  of 
the  shape  function  N.  is  that  its  value  is  unity  at  node  i  and 
is  zero  at  all  other  nodes.  For  the  quadrilateral  element,  the 
Nj  are  as  follows.  . 


"S^' 
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(2.27) 


^  =  1   (l+^)(l+Tl) 
■^    4 


N^  -  ^   (i-n(i+n) 

^    4 


^3  =  i  (i-n(i-n) 

■*    4 


^4  =  4  (1+^)(1-T1) 
^    4 


To  evaluate  the  stiffness  matrix  of  the  element,  we  need 
to  calculate  the  strain-displacement  transformation  matrix  for 
the  small  strain  case  as, 
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Yrz    = 
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As  the  interpolation  functions,  Nj ,  are  defined  in  the 
natural  coordinate  system  by  Equation  (2.27),  we  need  the 
relationship  between  r  and  z  derivatives  and  i  and  ^ 
derivatives.  By  using  chain  rules,  the  following  expressions 
can  be  obtained. 


31 


•  du 
dr 

du 
dz 


dw 

dr 

- 

.      = 

dw 

[   dz  \ 

dr 

dz 

dl 

dl 

dr 

dz 

dr\ 

^ 

dr 

dz 

di 

di 

dr 

dz 

dr\ 

dr] 

-1 

r  du  ] 

r  du  ] 

di 

^—1 

d^ 

. 

•  =  j^  ' 

du 

du 

[  ati  J 

[  an  J 

-1 

dw  ' 

dw 

a$ 

_._T 

dl 

■* 

dw 

'    =  J  ^  ^ 

dw 

[  an  J 

[  an  J 

(2.28) 


(2.29) 


In  the  above  equations,  J  is  the  Jacobian  operator 
relating  the  natural  coordinate  derivatives  to  the  local 
coordinate  derivatives.  The  inverse  of  Jacobian  matrix  must 
exist  in  order  for  the  isoparametric  element  to  be  used.  This 
means  that  there  is  a  one-to-one  correspondence  between  the 
natural  and  the  local  coordinates  of  the  elements.  In  actual 
computations,  the  sign  of  the  Jacobian  determinant  is 
monitored  at  the  Gaussian  integration  points  to  ensure  the 
existence  of  the  inverse  Jacobian.  Generally,  a  zero  or 
negative  value  of  the  Jacobian  determinant  indicates  input 
data  error  or  an  overly  distorted  element,  and  the 
computations  are  terminated  automatically  for  the  purpose  of 
node  renumbering  [22]. 

Equations  (2.28)  and  (2.29)  can  be  calculated  as  follows. 


dx 

dl 

dr 


=  1  (i+Ti)r,--i  (i+n)^2-|  (i-n)^3-^|  (i-n)^4 

=  1  il^^)  r,-l  (1-^)  r,-l  il-i)  r,-l  ilH)  r. 


(2.30) 
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9|  =  li^^r\)  z^-^[l^r\)  z^-^[l-r\)  z^^^{^-r\)  z, 

^  n^i)  z^^^  {i-D  z^-\  n-^)  z,-^  u^i)  z. 


dz 

di 

dz 

an 


(2.31) 


du 

as 

au 
an 


^  il+x])  u^-^  il^r])  u,-^  il-r^)  u,^^  il-^)  u, 
^  (1^^)  u,^^  (1-^)  u,-^  il-i)  u,-^  il^l)  u. 


(2.32) 


1^  =  1  il+x])  w^-^  {l+r\)w^-^  {l-ri)w^+^{l--q)w^ 


dw 
dw 


4 


dr)    4  4  4  ^ 


t:<sf-- 


(2.33) 


Eventually,  the  unknown  variables  of  the  system  will  be 
the  nodal  displacements  {  U  } 


{    U  }    =    {    u^    w^    U^    W^    U^    Wj    u^    w^  ]' 


(2.34) 


Then, 
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r  du  ] 

dr 

du 

4 

l+Ti    0    -d+ri)     0    -d-ri)     0        1-Ti        0 

1  +  5  0      1-5      0  -d-o  0  -d  +  o  0 

[  dz  J 

dw  ' 

dr 

dw 

4 

0    1+Ti    0    -d+r))     0    -d-Ti)    0        1-Ti 

_  0  1  +  $  0      1-$      0  -  d-o  0  -  d  +  o 

dz 

{  U} 


{  U} 


strain  can  be  expressed  with  nodal  variables  {  U  } ,  and 
also  stress  can  be  obtained  by  multiplying  the  strain  vector 
by  the  elasticity  matrix  [  C  ]g. 


{  e  }  =  [  B  ],  {  U} 


{   o  }    =    [   C]^    [  B]^   {   U) 


(2.35) 


In  the  axisymmetric  case,  unlike  other  isoparametric 
elements,  the  strain-displacement  matrix  [  B  ]^  has  terms 
containing  r.  Thus,  for  two  elements  of  identical  shapes  and 
material  properties  but  at  different  r  locations,  the  [  K  J^ 
will  be  different.  Also,  a  problem  arises  when  r  is  zero 
because  some  terms  have  r  in  the  denominator.  A  proper 
numerical  scheme,  such  as  Gaussian  quadrature,  can  avoid  this 
problem  [24] . 

In  the  axisymmetric  case,  {  e  }  becomes 
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{  e  }  =  {e^^  e^^  Yrz  ^ee^ 


(2.36) 


and  [  B  ]  becomes  4  by  8  matrix  as  follows, 


i+n 

0 

1  +  5 

(1  +  5)  (1+11) 

r 

0 

d+Ti) 

1+5 
0 

l+Tl 

1-5 

0 

(1-5)  d+Ti) 

r 

0 

d-Ti) 

1-5 

0 

-  d+n) 
-(1-5) 

0 

(1-5)  d+n) 

r 

0 

l-ri 

-(1-5) 

0 

-d-Ti) 
-(1  +  5) 

0 

(1  +  5)  d+Ti) 

r 

0 

-(1+5) 

i-n 

0 

(2.37) 


;■;*• 
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The  element  stiffness  matrix  can  be  formulated  by  using 
the  virtual  work  principle.  The  principle  requires  that  for  a 
body  under  equilibrium, 

f  be^'adV  =  j  ba^f^dV  +  j  ba^fsdS  ^  6a^f^  (2.38) 


where 


o   ;  actual  stress 

be    ;  virtual  strain  that  corresponds  to  the  imposed 

virtual  displacement 
6a  ;  virtual  displacement 
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fg  ;  actual  body  force 
fg  ;  actual  surface  force 
f^  ;  actual  nodal  force 


for  the  axisymmetric  element,  the  volume  integral  of  any 
function  A  becomes 


f  A  dv  =   f  A  rdrdzc^  (2.39) 


V. 


Volume  integration  can  be  written  in  the  natural 
coordinate  system.  Notice  that  r,  which  is  a  global 
coordinate,  appears  in  the  integration,  which  is  a  unique 
feature  of  the   axisymmetric  element. 


f  A  dv  =  2Tiff  A  r  det  [J]  d^dr\  (2.40) 


If  there  are  no  body  forces,   or  initial  stresses. 
Equation  (2.38)  becomes 

6a^|  [  B-\l[  C]^[  B],dV  {  U}   =  ba^  f^  (2. 41) 
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by  applying 


{  5e  }    =   [  B]^  {  bU) 
{   a  ]    =    [   C]^   {  e  }    =    [   C]^    [  B]^   {   U]  (2.42) 

[ba  }    =    [  N]     {bU] 


The   stiffness   matrix   is   obtained   as, 


[  i^],  =  /    [  B]l[  C]^[  B]^dV  (2.43) 


V. 


It  is  very  easy  and  convenient  to  use  Gauss-Legendre 
numerical  quadrature  integration  formula  to  evaluate  Equation 
(2.43).  The  Gauss-Legendre  integration  points  ^.  and  r|j  and 
corresponding  weighing  factor  a^^  are  found  in  Ref.  [24].  Care 
must  be  taken  to  substitute  r  in  Equations  (2.37)  and  (2.43). 
From  Equation  (2.26),  r  can  be  obtained  as 


=  ^  (l+$i)  (1+Ti^.)r,  +  -|  (1-?,)  {l+r]j)  r^ 

(2.44) 
■^1(1-$,)  (l-n^.)r3  +  -|(l  +  ^,)  (1-Ti,.)r, 


After  substituting  ^-  and  tIj  into  I    and  ti  at  the  previous 
equations,  [  K  ]  can  be  obtained  numerically  as 
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(2.45) 


The  assembly  of  element  stiffness  matrices,  applying 
boundary  conditions  and  the  solution  of  resulting  linear 
equations  follow  the  standard  procedures  that  can  be  found  in 
Ref .  [24,25] .  ; 

■■•^..-.-if 
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Figure  2.1  The  continuity  of  the  stresses  through  the 
thickness  of  the  composite  plate,  Ref.  [23] 
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Figure  2.2  The  continuity  of  the  strains  through  the 
thickness  of  the  composite  plate,  Ref.  [23] 
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Figure  2.3  Stacking  of  plies  to  one  element 
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Table  2.1  Material  Properties  of  AS4/3501-6 

(unit;  Gpa) 


El 

144.8 

^12=^13 

6.0 

^12=^13 

0.3 

^2=^3 


'23 


9.7 
3.6 

0.34 


*  Courtesy  of  Hercules  Advanced  Materials  and 
Systems  Company 


Axisymmetric  Element  Property  Matrix,  [  C  ]g 
for  [  0°,  ±45°,  90°  ] 

(unit  ;  Gpa) 


63.593    4.178 


4.178    11.105 


20.198    4.178 


4.500 


20.198 


4.17  8 


63.593 
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a  ;  Radius  of  the  plate 

c ;  Contact  radius 

R  ;  radius  of  the  indenter 

Gk  ;  Gap  between  indenter 


v,^:%>.^ 


Figure  2 . 4   Indentation  of  a  plate 
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Figure   2  .  5     Axisyitunetric  Element 


CHAPTER  3 
RESULTS  OF  FINITE  ELEMENT  ANALYSIS 

3 . 1  Introduction 
In  this  chapter,  the  accuracy  of  the  finite  element  -. 
program  is  evaluated  by  comparing  with  the  results  of  Sankar 
[10],  who  employed  a  finite  difference  method  to  solve  the  ; 
equilibrium  equations.  The  results  obtained  using  only  384 
finite  elements  (24  x  16)  in  the  present  study  compared  very 
well  with  that  of  Ref.[10].  In  Section  3.2,  available  methods  > 
of  increasing  toughness  and  damage  tolerance  are  briefly 
discussed.  In  Section  3.3,  the  finite  element  program  was 
applied  to  the  experimental  test  configuration  shown  in  Figure 
1.2  to  obtain  contact  force  vs.  center  deflection.  Without 
modification  of  the  total  stiffness  of  the  plate  due  to 
internal  damages,  finite  element  results  seem  to  deviate  from 
the  experimental  results.  These  deviations  can  be  attributed 
to  various  damages  as  will  be  discussed  in  Chapter  4.  However, 
in  the  elastic  range,  where  damage  is  minimal,  finite  element 
results  match  the  test  curves  well  (Figure  3.6).  In  Section 
3.4,   various  empirical   equations   including  the  Hertzian 
indentation  law  are  briefly  introduced  and  compared  with 
finite  element  results  in  the  elastic  range.  In  Section  3.4.1, 

44 
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the  finite  element  program  is  used  to  study  the  indentation  of 
solid  cylinders  and  plates  to  compare  with  Hertzian 
indentation  law.  Effects  of  plate  radius  and  indenter  radius 
on  indentation  and  deflection  of  the  plate  are  studied.  In 
Section  3.4.2,  Hertzian  contact  pressure  distribution,  which 
is  known  to  be  elliptic,  was  compared  with  finite  element 
h'."  results.  In  Section  3.4.3,  the  relation  between  contact  force 

•  V  and  contact  radius  was  studied.  It  is  known  that  in  a  half- 
space  contact  force  varies  as  the  cube  of  the  contact  radius. 
Finite  element  results  matched  very  well  with  the  available 
empirical  equations.  In  Section  3.4.4,  the  relation  between 
contact  force  and  center  deflection  was  studied.  In  Section 
3.5,  the  effect  of  friction  between  indenter  and  plate  was 
considered. 

3.2  Evaluation  of  the  Program 
Hybridization,  such  as  mixing  the  fiber  types,  thus 
utilizing  the  advantageous  properties  of  one  material  to 
overcome  deficiencies  of  another,  is  one  method  to  deal  with 
impact  damages.  The  higher  strain  to  failure  of  glass  fiber  or 
Aramid  fibers  (Kevlar)  led  to  substantial  increases  in 
threshold  energies  for  impact  damage,  with  the  energy  level 
required  to  produce  delamination  being  increased.  The 
Aramid/graphite/epoxy  hybridization  is  promising  because 
Aramid  material  is  relatively  inexpensive  and  it  increases 
impact  resistance  while  graphite  fibers  provide  compressive 
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strength,  which  Aramid  fibers  lack.  Tests  showed  that  doubling 
of  residual  strength  after  impact  was  possible  using 
Aramid/graphite  hybridization  compared  to  an  all  graphite 
specimen  [2]  . 

It  is  well  known  that  tough  resins  can  significantly    .  . 
reduce  the  damage  caused  by  impact  and  substantially  improve  **"  ^' 
the  residual  strength  following  impact.  This  is   because  the    ^  :i 
onset  of  delamination  and  the  delamination  growth  in  a     . 
composite  depends  on  its  interlaminar  fracture  toughness, 
which  in  turn  depends  on  the  toughness  of  the  matrix  material.     .     \ 
Widely  used  epoxy  resins   are  brittle,   but   crystalline 
thermoplastics  such  as  PEEK  (Poly-ether-ether-ketone)  can  be 
used  in  place  of  epoxy  to  increase  toughness.  Mode  I  fracture  ■ 
toughness  of  the  graphite/epoxy  ranges  80-200  J/m^  ,  while  the  '. 
fracture  toughness  of  graphite/PEEK  is  about  2500  J/m^  [2] 

Interleaving  is  a  process  in  which  thin  layers  of  ' 
thermoplastics  are  placed  in  between  layers  of  thermosetting 
prepreg  materials  to  form  a  tough  interlayer.  It  is  useful  if 
the  failure  occurs  through  delamination.  It  is  not  effective 
if  the  transverse  shear  is  high  and,  for  this  case,  an 
improvement  of  shear  modulus  of  the  resin  is  more  important 

Through  thickness  reinforcements,  such  as  a  three 
dimensional  braiding  process  or  through-the-thickness 
stitching,  can  improve  impact  resistance.  However,  in-plane 
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strength  will  be  reduced  because  of  wavy  fibers  and  irregular 
fiber  orientations. 

Sankar  [10]  performed  the  analysis  of  indentation  of 
hybrid  and  interleaved  plate  problems  by  using  a  finite 
difference  method  starting  from  the  axisymmetric  equilibrium 
equations. 


where 


(3.1) 


{  ^3  +  C,,  )  (  u,„  +  r-^u,^  )    +  C33iv^,^^ 
+  C44  {  w,^^  +  r'^w,^  )    =0 


(3.2) 
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^ee 
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(3.3) 


"^rz    ^4i    fzz 


(3.4) 


The  results  obtained  by  the  present  finite  element 
program  coincide  very  well  with  results  in  Ref.  [10].  Figure 
3.2  shows  contact  force  vs.  contact  radius  relationship  for 
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the  three  different  laminated  plates.  Material  properties  of 
the  composite  laminates  are  listed  in  Table  3.1. 

Figure  3 . 3  shows  interlaminar  shear  stress  in  the 
homogeneous  graphite/epoxy  laminate.  Figure  3.4  shows 
interlaminar  shear  stress  in  the  hybrid  laminate.  Employing 
the  finite  element  program  is  very  economical  because  only  384 
(24  X  16)  elements  were  used,  and  computing  time  was  about  30 
minutes  in  a  Vax  main  frame  computer.  Although  this  analysis 
does  not  account  for  the  internal  damages  of  the  plate,  it  . 
provides  information  about  the  possible  onset  of  delamination 
and  its  subsequent  propagation.  As  delamination  failure  occurs 
where  shear  stress  is  high,  shear  stress  distributions  such  as 
Figures  3.3  and  3.4  suggest  locations  where  delamination 
failures  possibly  may  initiate.  By  obtaining  other  strain  and 
stress  distributions  inside  plate,  possible  damage  pattern, 
such  as  matrix  crack  and  fiber  breakage,  and  its  extent  can  be  ■ 
predicted  conservatively  and  less  expensively  before  detailed 
damage  analysis.  Strengthening  by  inserting  tough  materials  or 
hybrid  laminates  to  minimize  delaminations  can  be  studied  by 
using  the  finite  element  program. 

Detail  damage  analyses  are  presented  in  Chapter  4. 


3.3  Comparison  with  Test  Results 
Jv.  The  finite  element  program  was  used  to  analyze  some  of 

ji.  •  ■-.        the  specimens  used  in  the  indentation  tests  discussed  in 
^  J-       Section  1.2.  The  contact  force  vs.  center  deflection  curve  was 
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obtained  without  considering  damages  inside  the  plate  and  was 
compared  with  the  test  results  in  Figure  3.5.  Even  in  the 
presumably  elastic  range,  finite  element  results  make  the 
plate  appear  stiffer.  However,  from  Hashin's  failure  criteria, 
which  will  be  discussed  in  Section  4.1,  the  damages  begin  to 
accumulate  at  about  140  N  contact  force  or  0.33  mm  center 
deflection.  Hence  the  actual  elastic  range  is  very  small. 
Comparing  the  results  only  in  this  elastic  range,  results  of 
finite  element  and  experimental  data  showed  fairly  good 
agreement.  A  total  of  512  elements  (32  x  16)  were  employed. 
Fine  mesh  was  used  near  the  center  of  the  plate,  and  coarser 
elements  were  employed  outside  the  contact  region.  Actually, 
finite  element  results  form  an  upper  bound  for  the  test 
results.  It  can  be  easily  explained  that  these  early 
deviations  of  the  stiffness  of  the  test  results  are  due  to  the 
inherent  imperfections  in  the  test  specimens  such  as  voids, 
debonds  and  nonuniform  fiber  orientations. 

3.4  Comparison  of  Contact  Behavior 
In  this  section,  various  results  concerning  indentation, 
contact  pressure  distribution,  contact  radius  and  center 
deflection  obtained  using  the  present  finite  element  method 
are  compared  with  available  results. 
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3.4.1  Hertzian  contact  law 

Because  of  the  extremely  difficult  nature  of  analytical 
solutions  of  contact  problem,  which  involves  permanent 
deformations  and  damage  such  as  matrix  cracks  and 
delamination,  an  experimental  approach  has  been  taken  to 
determine  the  contact  law  of  composite  materials.  Naturally, 
loading  and  unloading  curves  (contact  force-indentation 
relation)  are  different  because  damages  occur  even  at  low  load 
levels.  This  is  because  composite  materials,  unlike  monolithic  ■ 
materials,  have  very  low  resistance  to  the  damage  due  to  loads 
transverse  to  the  fiber  direction. 

For  the  loading  curve  of  contact  between  an  elastic 
sphere  and  elastic  half  space,  the  Hertzian  contact  law  [6]  is 


where  K^,  is  .  ^.; ,  *^ 


,^  .  1  ^,  :  l±_l^  .  1±_I^  yi  (3..) 


Subscript  s  denotes  the  indenter  and  t  denotes  the  target.  In 
Equation  (3.6),  R^  v^  and  E^  are  the  radius,  the  Poisson's 
ratio  and  the  Young's  modulus  of  the  indenter,  respectively. 
When  a  metal  indenter  like  a  steel  ball  indents  a  polymer 


-»■ . 
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composite,  the  indenter  can  be  assumed  rigid  because  Young's 
modulus  of  steel  is  much  higher  than  E^,  the  modulus  of 
elasticity  of  the  composite  in  the  transverse  direction.  For 
loading,  Sun  and  Yang  [6]  suggested  the  contact  force- 
indentation  relationship  in  the  form 


F  -  I  E^/r;  a3/2  (3,7) 


A  more  complex  form  of  contact  law  obtained  by  Conway 
[26]  for  an  elastic  spherical  impactor  striking  a  transversely 
isotropic  composite  medium  has  only  3  percent  difference  from 
Equation  (3.7)  in  the  example  problem  to  be  presented  in  this 
section.  In  Equation  (3.7),  C^^  of  [  C  ]  ^^  matrix  can  be 
replaced  by  E^  without  significant  error  (4  percent  difference 
in  this  example  problem) . 

The  unloading  formula  suggested  by  Sun  and  Yang  [6]  was 


F  =  -^  (a-ao)2-5   .  (3.8) 

a,       " 


where 


«o  =  %  [l-(^)  ^]   ,  a„>  a, 


a„  =  0  ,   a„  ^  a. 


(3.9) 


K.     .  ...  .  _.        "^ 


/ 
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In  Equation  (3.9),  a^  ^^  ^^^  indentation  corresponding  to  the 
maximum  contact  force  before  unloading.  It  is  assumed  that  the 
critical  indentation,  a^^.,  depends  only  on  the  material 
properties. 

Contact  force  vs.  indentation  graphs  were  obtained  by  the 
finite  element  program  for  the  clamped  solid  cylinder  of  20  mm 
radius  and  20  mm  thickness  of  a  hypothetical  isotropic 
material,  and  also  for  AS4/3501-6  composite  laminate  plate  of 
25.4  mm  radius.  The  material  properties  are  listed  in  Table 
3.1.  A  total  of  512  elements  (32  x  16)  were  employed.  Inside 
the  contact  region,  fine  mesh  is  used,  while  coarser  mesh  size 
is  assigned  outside  the  contact  region.  Each  computer  run  took 
75  minutes  on  the  average  in  the  Vax  main  frame  computer.  In 
the  case  of  indentation  of  the  thick  cylinder,  the  center 
displacement  at  the  top  becomes  indentation  directly.  For  the 
plate,  indentation  is  defined  as  displacement  difference 
between  top  and  bottom  surface  of  the  plate  at  the  center, 

a  =  w(0)    -  w{h)  (3.10) 

The  results  for  contact  force-indentation  are  presented 
in  Figures  3.8  and  3.9.  For  the  solid  cylinder  of  isotropic 
material,  the  finite  element  results  and  Hertzian  indentation 
law  of  Equation  (3.6)  agree  very  well  as  seen  in  Figure  3.8, 
although  errors  were  expected  because  we  have  not  solved  the 
half  space  contact  problem,  and  the  indenter  can  not  be 
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simulated  exactly  as  a  sphere  in  a  finite  element  program.  For 
the  plate  of  isotropic  material,  contact  forces  predicted  by 
the  finite  element  model  are  larger  than  forces  calculated  by 
Hertzian  contact  law  when  the  indentation  is  large  (Figure 
3.8).  Sankar  [27]  noticed  this  difference  in  the  contact 
coefficients  of  Equation  (3.5)  between  contacts  of  a  sphere 
with  a  half  space  and  with  a  plate.  He  suggested  that 
differences  of  the  definitions  of  indentation  between  half 
space  and  plate  may  be  the  cause  of  the  problem.  The  finite 
element  contact  forces  of  the  composite  plate  (4  mm  thickness) 
are  shown  in  Figure  3.9. 

To  check  the  influence  of  the  radius  of  the  plate,  2  5.4 
mm,  38.1  mm,  and  50.8  mm  radius  plates  are  compared  in  Figures 
3.10  and  3.11  (the  indenter  radius  is  12.7  mm).  Although  the 
contact  force  vs.  center  deflection  curve  shows  a  big 
difference  due  to  the  bending  effect  of  the  plate  (Figure 


3.10),  contact  force  vs.  indentation  curve  does  not  vary  so   -  [, 

much  (Figure  3.11).  At  small  indentation,  contact  force  seems 

to  be  independent  of  thickness  and  radius  of  the  target.        _  ;  ._ 

To  verify  that  contact  force  is  proportional  to  the 
square  root  of  indenter  radius,  two  different  indenter  radii, 
6.35  mm  and  12.7  mm,  were  used  to  indent  a  4  mm  composite 
plate  (Figure  3.12).  The  contact  force  curve  obtained  using 
the  finite  element  program  for  the  indenter  radius  6.3  5  mm  was 
denoted  as  F^  and  the  curve  for  12.7  mm  indenter  as  F^.    It  may 
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be  seen  that  for  the  small  indentation,  the  contact  force  is 
proportional  to  the  square  root  of  indenter  radius. 

3.4.2  Contact  pressure  distribution  of  the  plate 

The   contact   pressure   distribution   for   half-space 
indentation  is  as  follows  [27]. 


3  P 


a  =       ^   ^       Vl  -  {r/O^  (3.11) 

2  7t  C^ 


where  P  is  contact  force,  and  C  is  contact  radius.  Finite 
element  results  were  compared  with  these  equations  in  Figure 
3.13.  Obtaining  contact  stress  distribution  by  experimental 
test  is  not  easy  and  no  data  is  available  for  the  tests  done 
as  in  Section  1.2.  Contact  pressure  distribution  calculated 
by  the  finite  element  method  in  the  4  mm  plate  (Figure  3.13) 
is  elliptic.  Contact  pressure  of  the  finite  element  method  was 
computed  from  the  average  vertical  stress  o^  of  the  elements 
in  contact  with  the  indenter.  Although  calculated  contact 
forces  were  larger  than  predicted  by  the  half-space  formula. 
Equation  (3.11),  the  contact  pressure  distribution  is  still 
elliptical. 

3.4.3  Contact  radius  of  the  plate 

In  a  half-space  the  contact  force  is  proportional  to  the 
cube  of  the  contact  radius  C,  By  inserting  an  empirically 

..  ■       -  -^^  ■ 
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obtained  relationship  between  indentation  and  contact  radius 
for  the  indenter  radius  as  in  Equation  (3.12)  into  Hertzian 
contact  law,  Equation  (3.7),  we  can  obtain  Equation  (3.13) 
[28]. 


r2 
a  =  —  (3.12) 


P  =  ^   C^  (3.13) 

3^s 


Finite  element  results  for  a  4  mm  plate  are  compared  with  the 
above  equation  in  Figure  3.14.  In  Figure  3.15,  contact  force 
vs.  contact  radius  curves  of  plates  of  different  thicknesses 
are  compared.  As  plate  thickness  becomes  larger,  more  contact 
force  is  needed  to  maintain  the  same  radius  of  contact,  as 
expected.  However,  for  the  small  contact  radius,  contact  force 
can  be  assumed  to  be  independent  of  thickness  of  the  plate  as 
suggested  in  Equation  (3.13).         '   ;  ; 

3.4.4  Center  deflection  of  the  plate 

The  center  deflection  6  of  the  impacted  surface  of  the      A 
plate  can  be  expressed  as  the  sum  of  indentation  and  of      "' 
displacements  due  to  bending  and  shear.  In  the  case  of  small      • 
deflection,  membrane  effects  can  be  neglected.  Indentation  was 
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obtained  from  Equation  (3.7).  Plate  bending  stiffness,  Kp,  can 
be  expressed  as  [29], 


16-it  (1+v)  D 
(3  +v) a^ 


K^  =    ^-'^_^^-y>^^  (3.14) 


for  the  axisymmetric  plate  with  thickness  h  and  radius  a.  In 
Equation  (3.14),  D  is 


D  =   — (3.15) 

12  (l-v2) 


where,  E  can  be  replaced  by  C^^  of  the  matrix  [  C  ]^^  and  v, 
Poisson  ratio,  can  be  replaced  by  v^^.  Shear  deflection  can  be 
obtained  from 


T^^  =  G-^  =  — ^-  (3.16) 

''''  dr        2-Krh 


Plate  shear  stiffness,  K  ,  can  be  obtained  as 


K^  = ^ (3.17) 

^         i  Ln  a  -  Ln  C  ] 


and  C  can  be  obtained  from  Equation  (3.12). 

Then,  the  center  deflection  6  is  expressed  as 
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p\t-s^S..-^  (^-^'^ 


'-\±rK,^  K^ 


A     ™™  T-,1  at-e  are  compared  with 
^„4-  r-<=c;nlts  of  a  4  mm  piai:e  cii.-=     t- 
Finite  element  resuxT:i= 

.^atio.  <3..S,  in  n,.re  3..S.  Xt  is  .ncwn  ..at  B^at.on 
(3.»,  can  predict  ,00.  results  .or  the  elastic  response  o.     , 
..e  Plate.  This  result  su,,ests  t.at  center  ae.lection  .s  a 
„ore  reliable  parameter  t.an  indentation  to  obtain  ^^  <:--=!  .,,, 

force. 

*  .,  ■  —    <  /  ', 
T-S  Fri^^^•''^'T^  Effect  "'    'd. 

.    ^-         -i^  ci-nr5ied  bv  assuming  that  the 
The  effect  of  friction  is  studied  cy 

lar-rrfs   There  are  only 
^     -Fr-ir-t-ion  is  very  large.   iiit;j-'= 
coefficient  of   friction  ±  j  ^  . 

t.e  inaenter.  Stic,  condition  rec^ires  tnat  t.ere  .e  no 
relative  .otion  durin,  a  load  increment  between  indenter  and 
p,,te  at  t.e  contact  nodes.  T.is  condition  is  satisfied  .y 
constraining  t.e  radial  displacements  o.  t.e  nodes  a.ter  t.ev   . 
co.e  into  contact  wit.  t.e  indenter  t3C-31^  ..is  can  .e 
easily  implemented  in  t.e  finite  element  program  by  modify.n, 
t.e  contact  algorithm  w.ic.  was  developed  in  Section  . .  4 .  Now, 

.0  o9^  in  section  2.4  must  include 
displacements  in  Equation  (2.22)  m  Secti 

^^   r>-F  1-he  ut3per  nodes.  In  Equation  (2.22), 
radial  displacements  of  the  upper 

w  is  replaced  by 


{   w^}    =    i   Ui,w^i,U2'^2' "i-l'^i-l^ 


(3.19) 
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Treating  the  vertical  displacements  in  the  same  way  as  in 
Chapter  2.4,  we  now  put  zero  for  the  radial  displacements  in  : 
Equation  (3.19)  during  contact. 

Radial  force  is  very  large  compared  to  vertical  force 
because  of  the  high  modulus  of  radial  and  circumferential 
stiffness  due  to  the  presence  of  fibers.  Hence,  matrix  failure 
may  occur  virtually  immediately  after  contact  starts. 
Calculated  contact  force  vs.  center  deflection  for  both 
frictional  and  frictionless  contact  are  shown  in  Figure  3.17. 
The  effect  of  friction  is  significant  at  higher  contact 
forces.  However,  the  radial  stresses  developed  may  cause 
failure  even  at  smaller  loads,  and  the  result  shown  in  Figure 
3.17  may  not  be  applicable  in  practical  situations. 


■  ■i., 
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Rs  ;  Radius  of  indenter,  10  mm 

h  ;  Thickness  of  plate,  2  mm 

C  ;  Contact  radius 

a  ;  Radius  of  plate,  50  mm 


Figure  3.1  Plate  dimensions  in  numerical  examples 
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Table  3.1  Material  properties  of  the  plates  used  in  the 

examples  [10] , 

(unit,  GPa) 


Material 


E. 


E, 


1  50        7 

2  20        8 

3  3        3 


re 


3.8 

0.285 

0.265 

3.6 

0.270 

0.255 

1.1 

0.350 

0.350 

1.  Homogeneous  laminate  plate 

material  1,  thickness  2  mm 

2 .  Hybrid  laminate  plate 

material  1,  thickness  0.25  mm 
material  2,  thickness  0.2  5  mm 
material  1,  thickness   1.5  mm 

3.  Interleaf  laminate  plate 

material  1,  thickness  0.25  mm 
material  3,  thickness  0.25  mm 
material  1,  thickness   1.5  mm 
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Figure  3.2  Contact  force  vs.  contact  radius  relationship 
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Figure  3 . 3  Interlaminar  shear  stresses  in  the  graphite/epoxy 

laminates 
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Figure  3.4  Interlaminar  shear  stresses  in  the  hybrid 

laminates 
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Figure  3.5  Contact  force  vs.  center  deflection  curve 
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Figure  3.6  Contact  force  vs.  center  deflection  curve, 

elastic  range 
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Rs:  Radius  of  Indenter 
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Figure  3.7   Configurations  for  numerical  examples 
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Table  3.2  Material  properties  of  the  Numerical  applications 

Isotropic  Material 
E  =  2.0  GPa      v=0.25        G=   0.8  GPa 


[C],= 


E{l-v) 


(1+v) (l-2v) 


1-v 
0 

V 


1-v 
1 


0 

V 


1-v   1-v 


0 

l-2v 
2  (1-v) 

0 


1-v 

V 

1-v 
0 


[C],  = 


2.4  0.8    0  0.8 

0.8  2.4    0  0.8 

0  0  0.8  0 

0.8  0.8    0  2.4 


Composite  Material 


E     =   56.427   GPa        E,   =   10.688   GPa 
r  z 


^re   =    0.3 


G^^  =   4.5   GPa 


v^^  =   0.263 


[Ch 


63.593 

4.17  8 

0 

20.198 

4.17  8 

11.105 

0 

4.17  8 

0 

0 

4.5 

0 

20.198 

4.17  8 

0 

63.593 

J  V.'   i 


V  i. 
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Figure  3.8.   Contact  force  vs.  indentation  curve,  isotropic 

material 
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Figure  3.9.  Contact  force  vs.  indentation  curve,  composite 

plate  (4  mm  thickness) 
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Figure  3.10.  Effect  of  plate  radius  on  contact  force  vs, 

center  deflection  curve 
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Figure  3.11.  Effect  of  plate  radius  on  contact  force  vs. 

indentation  curve 


72 


c 
o 

o 

Z 
\^ 

o 
o 

o 

UL 

o 

(0 

c 
e 
o 


350 


300 


250 


200    - 


150    - 


100     - 


F1  (R  =  6.35  mm) 
F2  (R=  12.7  mm) 
/2"F1 


— I 1 1 — 

0.008  0.012 


"T" 


0.016  0.02 

Indentation  (mm) 


1 1 — 

0.024 


Figure  3.12.  Effect  of  indenter  radius  on  contact  force  vs, 

indentation  curve 
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Figure  3.13.  Contact  pressure  distribution 
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Figure  3.14.  Contact  force  vs.  contact  radius 
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Figure  3.15.  Contact  force  vs.  contact  radius,  comparison 
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Figure  3.16.  Contact  force  vs.  center  deflection 
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Figure  3.17.  Friction  effect 
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CHAPTER  4 
DAMAGE  ANALYSIS 

4.1  Failure  Criteria 

The  calculated  stresses  and  strains  from  the  axisymmetric 
elements  are  based  on  the  assumption  that  the  plate  is 
transversely  isotropic.  However,  the  plate  consists  of 
multiple  plies  where  each  ply  is  orthotropic  and  has  different 
fiber  orientation.  Hence,  direct  application  of  any  existing 
failure  criteria  which  are  based  on  stresses  in  the  material 
coordinates  is  not  possible.  An  estimate  of  stresses  in  the 
material  coordinates  has  to  be  obtained  from  the  axisymmetric 
results. 

The  existing  failure  criteria  for  composite  materials  can 
be  divided  into  two  types.  One  is  a  interactive  failure 
criterion  and  the  other  is  a  set  of  independent  failure 
criteria  for  different  failure  modes.  The  interactive  failure 
criterion  is  basically  an  empirical  technique,  and 
determination  of  the  components  in  the  polynomials  of  the 
criterion  requires  difficult  and  expensive  biaxial  tests.  The 
most  important  criterion  in  this  category  is  the  Tsai-Wu 
theory,  which  uses  a  general  quadratic  polynomial  including 
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linear  terms.  This  theory  says  that  a  failure  surface  in  the 
stress  space  exists  in  the  form  as 

F^Oj^  +  F^jO^Oj  =  1     ,     i,j=l,...,6        (4.1) 

where  F.  and  F..  are  second  and  fourth  order  strength  tensors. 
This  criterion  can  not  be  widely  used  because  of  difficulties 
in  obtaining  F...  These  criterion  only  predicts  the  onset  of 
failure  and  not  the  failure  modes.  However,  if  F,-  and  F^.j  can 
be  obtained  correctly,  this  criterion  may  be  satisfactory,  and 
it  is  simple  to  use  [13,14], 

Independent  failure  criteria,  typically  maximum  stress 
theory  or  maximum  strain  theory,  can  predict  the  onset  of 
failure  and  also  failure  modes.  Distinguishing  failure  modes 
is  very  important  for  the  design  of  composite  materials  and 
also  for  the  progressive  finite  element  program  which  must 
modify  its  total  stiffness  matrix  depending  on  the  failure 
modes.  However,  it  lacks  accounting  for  the  interaction  of 
various  components  of  stresses  or  strains  in  the  failure 
process.  ,   .^  i  j 

Good  failure  criteria  must  predict  the  principal  failure 
modes  separately,  such  as  fiber  breakages,  matrix  cracking  and 
delamination,  and  each  failure  mode  criterion  should  be 
interactive  such  as  the  quadratic  formula  while  they  do  not 
require  biaxial  tests  to  determine  the  coefficients.  The 
quadratic  failure  criteria  must  take  into  account  of  the 
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compressive  or  tensile  stresses.  The  work  of  Hashin,  which 
deals  with  the  three-dimensional  stresses  fits  this  category 
and  is  given  in  Ref.  [13].  It  leads  to  six  independent 
criteria,  given  in  Equations  (4.2)  to  (4.7).  As  input  to  the 
failure  criteria,  the  following  material  strength  parameters 
are  needed. 


X. 


'12 


'13 


'23 


X, 


fiber  direction  tensile  strength 

fiber  direction  compressive  strength 

transverse  tensile  strength 

transverse  compressive  strength 

shear  strength  in  the  1-2  plane 

shear  strength  in  the  1-3  plane 

shear  strength  in  the  2-3  plane 

interlaminar  (through-the-thickness)  tensile  stress 

interlaminar  shear  strength 


Fiber  failure  criteria  have  two  modes;  compressive  fiber 
buckling  mode  and  tensile  fiber  breakage  mode. 

For  tensile  fiber  direction  stress,  o^  >  0,  Hashin 's  fiber  ^,-f 
failure  criterion  is  ,     .  . 


X. 


'12 


\2 


^12  J 


'13 


^3 


1^7^ 


(4.2) 
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For  the  fiber  buckling  mode,  o^  <  0,  the  criterion  is 


\2 


X. 


=  1 


(4.3) 


c     I 


There  are  also  two  types  of  failure  modes  for  the  matrix 
cracking,  depending  upon  o^  +  O3,  such  that  when  (Og  +  O3)  >  0, 
the  matrix  cracking  criterion  is 


02+03 


)'-4' 


•23 


O2O3 


■12 


S^2 


■13 


I      ^13 


1     (4.4) 


and  when  (o,  +  O3)  <  0,  the  matrix  cracking  criterion  is 
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There  are  two  basic  approaches  to  predict  delamination 
initiation  [32].  One  is  the  mechanics  of  materials  approach, 
which  essentially  compares  the  local  state  of  stress  in  the 
interply  matrix  layer  where  delamination  occurs  with  relevant 
strength  parameters.  The  other  approach  is  based  on  the 
application  of  fracture  mechanics  principles  and  will  not  be 
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discussed  here.  However,  once  the  delamination  initiates,  the 
concept  of  a  critical  strain  energy  release  rate,  G^,  will  be 
applied  to  predict  delamination  growth. 

In  this  study,  delamination  was  predicted  by  the 
mechanics  of  materials  with  a  simple  quadratic  interactive 
formula  which  depends  upon  O3,  such  that  when  O3  <  0,  only 
shear  stresses  induce  delamination  such  as  Mode  II  type,  and 
the  criterion  is 
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when  O3  >  0,  it  behaves  like  a  combination  of  Mode  I  and 
Mode  II  type  crack  initiation,  and  the  criterion  is 


^]\2kAjk=i  .  (4.7) 


Here  X,  and  Sj  can  be  assumed  as  the  tensile  strength  and 
pure  shear  strength  of  the  matrix.  This  assumption  is 
justified  by  the  experimental  observations  that  toughened 
matrix  materials  apparently  reduce  the  onset  of  delaminations. 

Assume  damage  of  a  particular  nature  occurs  at  a  location 
given  by  a  in  the  0°  layer.  Then  the  same  type  and  extent  of 
damage  will  occur  in  the  45°  layer  at  (45°  +  a)  and  in  the  -45° 
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layer  at  (-45°  +  a)  and  so  on.  That  means  in  a  quasi-isotropic 
laminate  the  damage  pattern  will  look  like  a  stair-case  in  the 
form  of  a  helix.  This  phenomenon  has  been  observed  in  impact 
experiments  also  [15,33]. 

At  a  particular  segment  where  the  fiber  direction  is  at 
angle  6  to  the  x  axis,  the  stress  components  in  material 
coordinates  can  be  obtained  as  follows 
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(4.8) 


because  o^q    and  Og^  are  zero  for  the  axisymmetric  analysis, 
stresses  can  be  easily  obtained  as  in  Equation  (4.9). 


o^  =   cos^  6  Oj.  +   sin^  6  Oq 


O2  =  sin^  0  a^  +  cos^  0  Oq 


■12 


cos  6  sin  0  (  o_^  -  a^) 


(4.9) 


(^3  =  Ozz 


Ti3  =  COS  0  O^^ 


T23  =  -  sin  0  o^^ 


84 
By  combining  Equations  (4.2)  through  (4.7)  and  Equation 
(4.9),  the  failure  criteria  for  implementation  into  the  finite 
element  program  are  obtained  as  follows. 

Fiber  failure 

(i)         cos^e    o^   +      sin^e    Og      >      0 
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(    ii    )         cos^e    o^  +      sin^e    Og      <      0 


/  \2 

cos^6  Oj.  +  sin^6  Oq 


=  1 


Matrix   failure 
(i)  sin^e    o^   +   cos^e    Og  +   o^^      >      0 
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(ii) 


sin^e   o^  +   cos^e   Oq  +   o^^     <      0 
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Note  that  the  delamination  failure  criteria  are  not 
functions  of  6.  This  fact  can  also  be  verified  from  the  C-scan 
results  shown  in  Figure  1.5. 
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Figure   4 . 1    Partial    failure   patterns 
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Table  4.1.  Ultimate  strength  of  AS4/3501-6  lamina 

unit  (MPa) 


X^  ;  2172 
Y^  ;  53.8 
S,2  ;  120.7 


Sp-7  r     89.3 


S[   ;  89.3 


X(,   ;  1724 


Y(,   ;  2  05.5 


S^3  ;  120.7 


X, 


;  56.5 


*  Courtesy  of  Hercules  Aerospace  Company 
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4.2  Implementing  Delamination  into  the  Finite  Element  Program 
4.2.1  Introduction 

Delamination  initiates  from  where  the  shear  stress  is 
high,  such  as  above  the  mid-surface  of  the  plate  close  to  the 
contact  surface.  Then  delaminations  start  to  expand  radially, 
and  more  delaminations  appear  in  the  bottom  plies. 
Delaminations  can  be  simulated  in  the  finite  element  program 
by  creating  slip  surfaces  between  elements  as  seen  in  Figure 
4.2.  New  nodes  are  created,  and  rigid  elements  (spring 
elements  with  very  large  spring  constants  K)  are  added  between 
two  nodes  across  slip  surfaces  when  the  normal  force  is 
compressive  to  prevent  interpenetration  of  the  nodes.  Adding 
new  surfaces  and  consequently  adding  new  nodes  makes  a 
cumbersome  procedure,  which  is  a  weakness  of  the  finite 
element  method.  The  total  stiffness  matrix  of  the  plate  must 
be  recalculated  after  delaminations  occur  or  grow.  Also  the 
dimension  of  the  total  matrix  is  increased  due  to  the 
introduction  of  new  nodes.  If  tensile  stress  exists  at  the  tip 
of  delaminations,  the  spring  element  must  be  eliminated  to 
allow  Mode  I  type  (Opening  Mode)  delamination  propagation. 
Then,  iteration  of  the  program  is  necessary  to  assure  that  the 
rigid  elements  are  placed  only  where  compressive  stresses 
exist.  However,  unless  the  contact  force  is  high, 
delaminations  seemed  to  appear  only  at  the  upper  portion  of 
the  plate  where  through-the-thickness  stresses  o^  are 
compressive.  The  present  finite  element  program  monitors  the 
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signs  of  all  o  across  the  delamination  surfaces  while 
implementing  Hashin's  failure  criteria,  and  it  was  found  that 
all  o  's  are  compressive  in  the  present  study.  It  is  safe  to 
assume  that  in  this  study  the  delamination  growth  is  Mode  II 
type  (Shear  Mode)  and  to  put  rigid  elements  at  all  the  nodes 
across  delamination  surfaces.  When  delaminations  become 
extensive,  especially  beyond  the  mid-surface  of  the  plate, 
then  the  tip  of  the  delaminations  could  have  high  tensile 
stresses.  This  Mode  I  crack  growth  combined  with  Mode  II  type 
will  accelerate  delamination  growth  faster  than  by  Mode  II 
delamination  alone. 

4.2.2  Numerical  Experiment 

A  simple  numerical  experiment  was  done  to  check  the 
effect  of  delaminations  on  the  contact  force  vs.  center 
deflection  curve  and  also  to  see  the  frictional  effect  across 
the  delamination  surface.  Assume  that  the  experimental  test 
configuration  shown  in  Figure  1.2  has  pre-existing 
delaminations  of  12.7  mm  diameter  uniformly  distributed 
through  the  thickness.  As  16  elements  were  employed  through 
the  thickness,  15  delaminations  were  assumed  in  the  finite 
element  program.  An  incremental  point  force  was  applied  at  the 
center  of  the  plate.  Rigid  elements  were  placed  across  all  the 
delamination  surface  nodes  as  shown  in  Figure  4.3. 

Figure  4.4  shows  contact  force  vs.  center  deflection 
curve  of  the  plate  with  pre-existing  delaminations  compared 
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with  experimental  test  results.  Delaitiinations  reduced  the 
slope  of  the  contact  force  vs.  center  deflection  curve.  It  is 
apparent  that  this  numerical  experiment  has  exaggerated  the 
extent  of  delaminations.  However,  it  justifies  the  conclusion 
that  delamination  failures  can  change  the  contact  force  vs. 
center  deflection  curve.  Combined  with  matrix  failures, 
delamination  failures  will  change  the  slope  of  the  contact  vs. 
center  deflection  curve  shown  in  Figure  3 , 5  which  was  obtained 
for  the  elastic  range  by  the  finite  element  program. 

4.3  Implementing  Fiber  and  Matrix  Failure  into  the 
Finite  Element  Program 

When  the  matrix  fails,  the  elastic  modulus  of  the  matrix, 
E  .  can  be  replaced  by  zero,  so  that,  from  elementary 
composite  material  theory, 


^1  =  EfV,  +   £„v,  =  EfVf 


Because  the  value  of  E^v^j^  is  very  small  compared  with  E^v^,  the 
original  E^  can  be  used  after  matrix  failure.  However,  the 
other  elastic  moduli,  E^    and  E3,  and  Poisson's  ratios  become   ,. 
zero.  With  these  new  values  input,  [  C  j^^  computed  by  the 
method  described  in  Section  2.3  becomes  singular  with  zeros  ^^  r  i  mJ 
for  C-,^  and  C,, .  Also  the  total  stiffness  matrix  assembled  with  .  ^^,  ^ 
these   singular   element   stiffness  matrices   can   also  be 
singular.  It  is  found  that  the  singularity  can  better  be 
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avoided  by  assigning  very  small  values  to  the  zero  diagonal 
terms  in  the  [  C  J^^  matrix  rather  than  to  the  diagonal  terms 
of  the  total  stiffness  matrix. 

Matrix  failures  were  initiated  first  and  then  expanded 
extensively.  There  were  no  fiber  failures  for  the  range  of 
calculations  in  the  present  study. 

The  finite  element  program  could  easily  handle  matrix 
failures  by  assigning  reduced  [  C  ]g^  values  corresponding  to 
damaged  elements  which  were  detected  by  Hashin's  failure 
criteria.  When  partial  failures  are  detected,  [  C  ]  ^^  must  be 
reduced  proportionally  depending  upon  the  degree  of  the 
failure.  However,  in  this  research,  partial  failures  were 
detected  only  at  a  small  portion  of  the  fringes  of  totally 
failed  elements.  So  partial  failures  could  be  either  ignored 
or  considered  as  total  failure  depending  upon  the  degree  of 
partiality  without  losing  much  accuracy.  The  total  stiffness 
matrix  is  reassembled  and  the  contact  algorithm  is  applied 
again  according  to  the  method  discussed  in  the  next  section. 

4.4  Incremental  Finite  Element  Method  for  Damage  Analysis 
Hashin's  failure  criteria  discussed  in  Section  4.1  have 
been  implemented  in  the  finite  element  program.  After  all  the 
displacements  are  calculated  as  in  Equation  (2.18),  stresses 
in  each  element  can  be  calculated  according  to  Equation 
(2.35).   The  finite  element  program  checks  Hashin's  failure 
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criteria  with  the  element  stresses  and  prints  damage 
information. 

When  damages  are  detected,  the  total  stiffness  matrix  is 
reassembled  to  take  account  of  delaminations  and  matrix 
failures  as  discussed  in  Sections  4.2  and  4.3.  When  failures 
occur  inside  the  plate,  the  displacements  of  all  nodes  except 
the  contacted  nodes  beneath  the  indenter  are  changed.  Among 
the  accumulated  displacements,  only  the  vertical  displacements 
of  the  contacted  nodes  of  the  plate,  d^,  remain  unchanged  and 
become  the  input  to  Equation  (4.10)  which  solves  the  mixed 
boundary  value  problem.  Here,  d^.  becomes  the  prescribed 
displacement  boundary  condition,  and  zero  traction  forces  at 
the  uncontacted  nodes  become  zero  force  boundary  conditions. 
Other  displacements  d^,  and  new  contact  forces  F^  can  be 
obtained  from  Equations  (4.11)  and  (4.12)  respectively. 
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The  incremental  contact  algorithm  is  then  applied  to  close  the 
gap  at  the  next  candidate  contact  node  as  discussed  in  Section 
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2,4  until  new  damages  are  developed  or  until  the  assigned 
contact  radius  is  reached. 

Sometimes,  stiffness  reduction  due  to  damage  causes  a 
node  just  outside  the  contact  region  to  penetrate  into  the 
indenter.  Since  this  is  not  physically  possible,  the  program 
imposes  a  displacement  boundary  condition  on  such  nodes  so 
that  they  just  touch  the  indenter,  and  the  problem  is  solved 
again  using  the  new  stiffness  matrix.  Equation  (2.25)  must  be 
modified  as 


W  ■    —    W  ■    +    Q- 
^     w^ w^ y^  (4.13) 

^     (  1  -d^  ) 


to  maintain  smooth  contact  as  seen  in  Figure  4.5.  When  6,.  is 
negative,  the  indenter  must  touch  that  node  automatically  and 
proceed  to  the  next  contact  node  to  maintain  smooth  contact. 
Here,  Gap  g.  is  the  z-distance  of  the  indenter  from  the 
corresponding  two  contact  candidate  nodes,  as  shown  in  Figure 
4.5. 

Usually,  matrix  failure  occurs  beneath  the  indenter 
because  of  compression.  At  almost  the  same  time,  the  matrix  at 
the  bottom  of  the  plate  starts  to  fail  because  of  the  bending 
of  the  plate.  Delamination  occurs  near  the  top  surface  of  the 
plate  beneath  the  indenter  because  of  high  shear  stresses,  as 
seen  in  Figures  3.3  and  3.4,  and  in  Figure  4.8. 
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When  matrix  cracking  and  delamination  failure  occur,  the 
slope  of  the  contact  force  vs.  center  deflection  curve  changes 
smoothly,  as  seen  in  experimental  data  (Figure  3.5).  It  is 
speculated  that  the  load  cell  attached  at  the  indenter  may  not 
catch  the  slight  load  drop  due  to  matrix  failure,  which 
behaves  like  plastic  deformation  [9].  The  finite  element 
program  recalculates  the  new  contact  force  and  center 
deflection  with  updated  damage  information  and  continuously 
interpolates  with  the  point  of  the  previous  damage  status. 
Although  contact  force  drops  with  increased  internal  damages, 
as  shown  by  the  vertical  dotted  lines  in  Figure  4.6,  it  is 
reasonable  to  draw  the  contact  force  curve  continuously  with 
respect  to  center  deflection  as  a  solid  line,  as  in  Figure 
3.5.  Contact  force  was  assumed  to  proceed  ignoring  the  effects 
of  small  growth  of  internal  damage  in  Figure  4.6.  This  does 
not  explain  the  sudden  drop  of  contact  force  shown  in  Figure 
1.3.  In  Ref.  [9],  this  sudden  drop  was  assumed  to  be  largely 
due  to  fiber  failure.  However,  examination  of  specimens  of  the 
impacted  or  indented  plates  at  the  University  of  Florida  did 
not  reveal  fiber  failures  [11].  Damage  patterns  of  the 
composite  plate  obtained  from  the  finite  element  program  are 
shown  in  Figures  4.8-4.10,  where  dotted  area  represents  matrix 
failed  elements  and  delaminations  are  shown  as  solid  lines. 
Judging  from  the  failure  growth  map  (Figures  4.8-4.10) ,  it  can 
be  assumed  that  extensive  matrix  failures,  which  drastically 
reduce  the  transverse  stiffness,  combined  with  delamination 
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growth,  induce  this  sudden  drop  of  contact  force. 
Delaminations  are  expected  to  initiate  from  matrix  cracks,  and 
including  this  initiation  phenomenon  in  the  finite  element 
program  requires  further  research. 

By  changing  element  mesh  sizes,  three  contact  force  vs. 
center  deflection  curves  were  obtained,  as  shown  in  Figure 
4.7,  where  they  are  compared  with  the  experimental  curves 
reproduced  from  Figure  3.5.  In  Figure  4.7,  the  curve  C  has  the 
coarsest  mesh  size,  and  the  curve  A  which  is  the  reproduction 
from  Figure  4.6  has  the  finest  mesh  size. 

Although  the  finite  element  program  appeared  to  be  an 
upper  bound  of  the  test  results  in  the  elastic  range,  it  was 
no  longer  an  upper  bound  when  the  damage  analysis  was 
included.  Also,  the  damage  analysis  was  highly  path  dependent. 
Change  of  mesh  size  may  alter  output,  and  results  can  vary  as 
much  as  the  experimental  data.  The  finite  element  program  was 
stopped  at  near  2000  N  contact  force  or  0.7  mm  center 
deflection  because  of  two  reasons.  First,  some  of  the  top 
layer  elements  of  the  plate  which  the  indenter  had  not 
contacted  yet  were  damaged.  Recalculated  contact  force  became 
near  zero  because  the  damaged  top  surface  did  not  give 
resistance  to  indentation.  Damage  on  the  top  surface  was  very 
sensitive  to  the  contact  force.  Second,  the  stiffness  matrix 
of  the  system  became  singular  numerically.  Assumption  of 
r  C  1   matrix  for  the  element  in  which  matrix  has  failed  as 
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discussed  in  Section  4.3  must  be  too  conservative.   By 
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increasing  values  of  C^^  and  C33,  the  numerical  singularity  can 
be  avoided.  By  numerical  experiments,  good  values,  which  can 
make  finite  element  results  match  well  with  experimental  data, 
might  possibly  be  found. 
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Figure  4.2  Delamination  simulation 
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Figure  4.3  Numerical  application  of  delamination  effect 
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Figure  4.4  Delamination  effect  on  Contact  force  vs.  center 
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Figure  4.5  Geometry  of  contact  surface 
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Figure  4 . 6  Contact  force  vs .  center  deflection  curve 
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Figure  4.7  Contact  force  vs.  center  deflection  curve 
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Figure  4.8  Damage  pattern  of  the  composite  plate 
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Figure  4.9  Damage  pattern  of  the  composite  plate 
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Figure  4.10  Damage  pattern  of  the  composite  plate 


CHAPTER  5 
CONCLUSIONS 

5.1  Conclusion 
The  contact  problem  of  a  circular  composite  plate  by  a 
spherical  indenter  was  analyzed  by  a  finite  element  program 
which  employed  quadrilateral,  isoparametric,  axisymmetric 
elements.  After  averaging  the  orthotropic  laminar  material 
properties,  axisymmetric  element  properties  of  cross-ply  or 
quasi-isotropic  laminate  turned  out  to  be  transversely 
isotropic.  Hashin's  failure  criteria,  which  were  written  in 
material  coordinates,  were  transformed  to  be  applicable  to  the 
axisymmetric  element.  Delamination  criteria  turned  out  to  be 
independent  of  angle  transformation.  A  circular  pattern  of 
delamination  was  also  shown  from  the  C-scan  picture  of 
experimental  data.  For  the  AS4/3501-6  graphite/ epoxy  plate 
(25.4  mm  radius),  the  elastic  region  is  only  up  to  140  Newton 
contact  force  or  0.3  3  mm  center  deflection  for  a  spherical 
indenter  of  radius  12.7  mm.  In  this  elastic  region,  results  of 
the  finite  element  program  matched  well  with  the  experimental 
data  in  most  cases,  and  were  an  upper  bound  of  the  test  data 
for  the  remaining.  The  finite  element  analysis  showed  that 
matrix  failure  was  initiated  both  from  top  and  bottom  - 
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beneath  the  indenter  because  of  compression  and  in  the  bottom 
surface  of  the  plate  due  to  bending.  Matrix  failure  expanded 
radially  and  also  toward  the  mid-surface  of  the  plate  from 
both  the  top  and  bottom.  Delamination  initiated  at  almost  the 
same  time  as  matrix  failure  near  the  top  surface  beneath  the 
indenter.  Delamination  expanded  radially  and  also  downward. 
When  the  delamination  and  matrix  failures  were  accounted  for, 
the  finite  element  results  matched  well  with  experimental  data 
in  predicting  the  force-deflection  relationship,  as  shown  in 
Figure  4.7.  Combined  extensive  matrix  failure  and 
delamination,  which  make  it  very  compliant  at  the  center  of 
the  plate,  is  believed  to  be  the  cause  of  the  sudden  drop  of 
the  contact  force  at  the  contact  force  vs.  center  deflection 
curve.  Widely  used  contact  laws  and  equations  are  compared 
with  finite  element  results. 

5 . 2  Comments  for  Future  Research 
The  sudden  drop  of  the  contact  force  could  not  be 
explained  without  modeling  the  delamination  propagation,  which 
is  very  complex.  Connecting  nodes  with  vertical  and  horizontal 
springs  and  iterating  schemes  to  account  for  the  growth  of 
delamination  are  needed.  Also,  as  matrix  cracks  will  lead  to 
delaminations,  understanding  of  this  growth  is  needed.  When 
the  matrix  fails,  a  better  model  of  the  constitutive  equations 
of  the  finite  element  is  required  for  accurate  prediction  of 
further  damage  propagation. 
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APPENDIX 

Description  of  the  Program 
Because  of  a  computing  time  limitation,  this  program 
assembles  and  condenses  the  total  stiffness  matrix  one  time 
i':'        only.  A  totally  automatic  program  for  a  damage  analysis  can  be 
^\  - •      developed  by  modifying  inputs  and  outputs  of  this  program. 
1.  Read  input  data 

A.  Total  contact  candidate  node  number  (NCAND) 

B.  Total  number  of  nodes  (NODES) 

Total  number  of  elements  (NELMS)  ■< 

Number  of  different  materials  (MPTOTAL) ;  reduced 
1  modulus  matrices  due  to  matrix  and  fiber  failure 

: ":  can  be  applied  here  j 

C.  Number  of  total  radial  elements  (NA) 

Number  of  radial  fine  mesh  elements  (NAFINE)  ..    ■'  ■ '.M- 

Number  of  elements  through  the  thickness  (NH)  '  '; 

Radius  of  the  plate  (AR)  ~        ;..^ 

Thickness  of  the  plate  (HT)  '! 
Radius  of  total  fine  mesh  elements  (ARFINE) 

D.  Boundary  conditions  (see  program  list)  •  '  A^ 

E.  Material  identification  number  (MP)  '\ 
'             F.  Material  properties  (Cll,  C12,  C13,  C33,  C44) 

G.  Ultimate  strength  of  the  materials  (XT,  XC,  XI,  YT, 
YC,  SI,  S12,  S13,  S23) 

H.  Indenter  radius  (RB) 

I.  Maximum  node  identification  number  to  calculate  ,  •>.   ' 
Contact  process  (NTRY) ;  when  indenter  contacts 
node  NCONT(NTRY) ,  program  stops.  ^^^ 

J.  Restarting  identification  number  (NREST)-  last  vi 

contacted  node  identification  number  to  *-^ 

recalculate  contact  force  due  to  change  of  •  /.4^' 

internal  damage;  zero  for  elastic  analysis        .^,     **' 

K.  Total  number  of  delamination  layers  (NTCR) 

L.  Identification  of  each  delamination  layer 

(NCRACK(I));  number  of  delaminated  elements        >     '>' 
, - ;  (MECR(J,I);  first  element  identification 

5        ;,,   M.  Total  number  of  matrix-failed  layers  (NRMF) 
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NRMFl;  number  of  elements  at  each  layer 

NRMF2 ;  first  node  identification 
N.  Total  number  of  fiber-failed  layers  (NRFF) 

NRFFl;  number  of  elements  at  each  layer 

NRFF2 ;  first  node  identification 
O.  Vertical  displacements  of  contacted  nodes,  input  for 

the  damage  analysis;  (DISP{I) , I=1,NREST) ) 

2.  Automatically  generate  elements,  nodes,  element  material 
properties  information. 

3.  If  restart  identification  number  is  not  zero  (step  1,  J), 
then 

A.  Create  new  delamination  surface  by  generating  new 
nodes,  solid  element  across  free  surface  to  prevent 
interpenetration  of  the  nodes. 

B.  Assign  new  reduced  stiffness  matrix  to  the  failed 
(matrix,  fiber)  elements. 

4.  Construct  global  stiffness  matrix 

5.  Create  condensed  stiffness  matrix  as  in  Equation  (2.18) 

6.  If  restart,  solve  Equation  (4.10),  then  proceed  to  next        ..■. 
node  contact  -^^  j^l 

7.  Create  sub-condensed  stiffness  matrix  as  Equation  (2.22) 

8.  Apply  contact  algorithm  as  in  Equation  (2.25)  or  Equation 
(4.13) 

9 .  Obtain  displacements  and  forces 

10.  Obtain  strains  and  stresses 

11.  Apply  modified  Hashin's  failure  criteria 

12.  Print  failure  information 

13.  Repeat  step  7  -  step  12  until  assigned  node  (step  1,  I) 
is  contacted. 


List  of  the  Finite  Element  Program 

PROGRAM  CONTACT 

IMPLICIT  REAL*8(A-H,0-Z) 
C    INPUT  INFORMATION 
C    1)   NCAND  /  15 
C    2)   NA,  NAFINE,  NH,  AR,  ARFINE,  HT  /  314,  3F6.3 
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C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 


3) 
4) 


5) 


6) 
7) 


8) 

9) 
10) 

11) 
12) 
13) 
14) 
15) 
16) 
17) 


/3I5 

OF  FIXED 


D.O.F. 


RADIUSWISE 


NELMS,  NODES,  MPTOTAL  /  315 
BOUNDARY  CONDITIONS  /  N1,N2,N3 

Nl;  FIRST  NODE  IDENTIFICATION 

N2;  DIRECTION,  R=l,  Z=2 

N3;  INCREMENTAL  NUMBER 

POSITIVE  :  THICKNESSWISE 
NEGATIVE  :  FIX  ALONG  THE 

END  OF  INPUT;  PUT  0  IN  Nl 
ASSIGNING  MATERIAL  PROPERTY  /  N1,N2,N3  /  315 

Nl,  FIRST  ELEMENT  IDENTIFICATION  OF  THE  LAMINAR 

N2,  MATERIAL  PROPERTY  IDENTIFICATION  OF  THE  LAYER 

N3,  NUMBER  OF  STACKS  OF  SAME  MATERIAL  LAMINATES 

END  OF  INPUT;  PUT  0  IN  Nl 
MP,  Cll,  C12,  C13,  C33,  C44  /  15,  5(F12.4) 


{S, 


y  = 


[  C,..  ]  {  E,  } 


REAd"^  ULTIMATE  STRENGTH  OF  i-HE  MAI^ERIAL  /  3F10.3 

XT,  XC,  XI 

YT,  YC,  SI 

S12,  S13,  S23 
RB  /  F10.5 
NTRY  /  15 
NREST  /  15 
NTCR  /  15 
NCRACK(I) 
NRMF  /I5 

NRMFl,  NRMF2  /  15, 
NRFF  /  15 

NRFFl,  NRFF2  /  15,  15 
DISP(I)  /  F20.12;  (1=1, NREST) 
OUTPUT  (DELETED  IN  THIS  LIST) 

MATRIX  AND  FIBER  FAILURE  (KFAIL(I)) 


MECR(J,I)  /  15,  15 


15 


-'/» 

: 


1  ;  TOTAL  MATRIX  FAILURE 

2  ;  PARTIAL  MATRIX  FAILURE 

3  ;  TOTAL  FIBER  FAILURE 

4  ;  PARTIAL  FIBER  FAILURE 
DELAMINATION  (KDELA(I) ) 

1  ;  MODE  II  FAILURE 

2  ;  MODE  I  AND  MODE  II,  TENSILE  STRESS  (Sz)  EXITS 
DIMENSION  RC(650) ,ZC(650) ,NBC(650,2) ,MNODE (550, 4) 
DIMENSION  S(8,8) ,CORD(2,4) ,D(4,4) ,B(4,8) ,DB(4) ,H(4) 
DIMENSION  GAUSS(4) ,WEIGHT(4) ,P(2,4) ,XJI(2,2) ,XJ(2,2) 
DIMENSION  MPROP(550) ,C11(5) ,C12(5) ,C13(5) ,C33(5) 
DIMENSION  C44 (5) ,NCONT(20) ,AC(20,20) ,STREC(4,2) 
DIMENSION  AX(20,20) ,BX(20,1250) , CX (1250 , 1250) 
DIMENSION  AA(20,20) ,AB(20,20) ,ZM(20) ,FXX(20) 
DIMENSION  WKAR(5500) ,RF(20) ,ZR(1250) ,ZU(20) ,ZZ(1,1) 
DIMENSION  DISP(1250) ,F(20) ,X(1250,20) ,MECR(30,15) 
DIMENSION  QUE (8) ,NSOLID(150, 2) ,ST(4) ,NCRACK(15) 
DIMENSION  STRESS(550,4) , STRAIN (550 , 4) ,SA(4) ,SN(4) 
DIMENSION  KFAIL(550),  KDELA (550) , ASN (4) 

DATA  GAUSS/0. 5773502692D0,-0.5773502692D0,0. ,0./ 
OPEN (UNIT=1 , FILE= ' ZI . DAT ' , STATUS= ' OLD • ) 
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OPEN (UNIT=2 , FILE= ' ZO . DAT ' , STATUS= ' NEW ' ) 
C   ****  READ  NODE  AND  ELEMENT  DATA 
READ (1,2000)  NCAND 
2000  F0RMAT(I5) 

WRITE (2, 1000)  NCAND 
1000  FORMAT (1 OX, 'NUMBER  OF  CONTACT  CANDIDATES ', 5X, 14) 

READ(1,2010)  NA,NAFINE,NH,AR,ARFINE,Ht 
2010  FORMAT(3I4,3F6.3) 

WRITE (2,1010)  NA , NAFINE , NH , AR , ARFINE , HT 
1010  FORMAT(/, 3 X, 'NUMBER  OF  DIVISION  OF  TOTAL  RADIUS ', 14 ,/ , 
+ 3 X, 'NUMBER  OF  FINE  DIVISION' , 14 ,/, 3X, ' DIVISION  OF 
+THICKNESS' ,I4,/,3X, 'TOTAL  RADIUS' ,F10.5,/,3X, 'FINE 
+MESHED  RADIUS' ,F10.5,/,3X, 'THICKNESS' ,F10.5,//) 
IC=0 

READ (1,2020)  NELMS , NODES , MPTOTAL 
2020  FORMAT(3I5) 
C  ****  AUTOMATIC  NODE  GENERATION 
DO  1  N=l, NODES 
DO  1  1=1,2 
NBC(N,I)=1 
1   CONTINUE 
C 

YH=NH 

YFINE=NAFINE 
XAFINE=ARFINE/YFINE 
YAC=NA-NAFINE 
XAC= ( AR-ARFINE ) / YAC 
XH=HT/YH 

DO  2  N=1,NAFINE+1 
I=N-1 
XI=I 

NI=(NH+1)*I+1 
RC(NI)=XAFINE*XI 
ZC(NI)=0.0 
IAUT0=1 
100  RC(NI+IAUTO)=RC(NI) 

ZC (NI+IAUTO) =ZC (NI+IAUTO-1) +XH 
IAUT0=IAUT0+1 

IF(IAUT0.NE.NH+1)  GO  TO  100 
2  CONTINUE 
C 

DO  4  N=1,NA-NAFINE 
I=N-1 
XI=I 

NI=(NH+1)*I+1 
NQX= (NH+1) * (NAFINE+1) 
RC (NI+NQX) =XAC* (XI+1 . ) +ARFINE 
ZC(NI+NQX)=0.0 
IAUT0=1 
104  RC(NI+NQX+IAUTO)=RC (NI+NQX) 

ZC(NI+NQX+IAUT0)=ZC(NI+NQX+IAUT0-1)+XH 
IAUT0=IAUT0+1 


KSfyr^ -  ■-r.y^'^r^       --  ■_  ■  ■■   ■■  ■::  jr' -  :-  ■  '%>^j^' 
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IF(IAUT0.NE.NH+1)    GO   TO    104 

4  CONTINUE 
C  ****  READ  FIXED  DEGREE  OF  FREEDOM 

K.         C       N;FIRST  NODE  NUMBER  TO  BE  FIXED 
:'         C       I;  DEGREE  OF  FREEDOM  TO  BE  FIXED  (  R=1,Z=2  ) 
V        C       L;  INCREMENTAL  NUMBER.  0  FOR  SINGLE  NODE  INPUT 
C       PUT  N=0  AT  THE  END  OF  THE  DATA 
2030  FORMAT(3I5) 
fl'  101  READ  (1,2030)  N,I,L 

■   .  IF(N.EQ.O)  GO  TO  102 

NBC(N,I)=0 

IF(L.EQ.O)  GO  TO  101 
x'"-  IF(L.LT.O)  GO  TO  103 

DO  3  J=1,L 
«    ^  3  NBC(N+J,I)=0 

GO  TO  101 
[-  103  DO  J=1,NA 

^  NI=(NH+1)*J+N 

NBC(NI,I)=0 
■ ■  END  DO 

GO  TO  101 
102  WRITE(2,1020) 
1020  FORMAT (/, 3 X, 'AUTOMATIC  GRID  AND  ELEMENT  GENERATION',/) 
C  ****  AUTOMATIC  ELEMENT  GENERATION 
•  MN0DE(l,l)=NH+3 

MNODE(l,2)=2  "j^ 

MN0DE(1,3)=1  '*' 

MNODE (1,4) =NH+2 

DO  5  N=1,NA 
'.  .  I=N-1 

XI=I 

NI=NH*I+1 

DO  6  J=l,4 
' ">  MNODE (NI , J) =MNODE ( 1 , J) + (NH+1) *I 

6  CONTINUE 
IAUT0=1 

110  DO  7  L=l,4 

MNODE (NI+IAUTO, L) =MNODE (NI+IAUTO-1 , L) +1 

7  CONTINUE 
IAUT0=IAUT0+1 
IF(IAUTO.NE.NH)  GO  TO  110 

5  CONTINUE 
'r-:                    c    ****    AUTOMATIC  MATERIAL  PROPERTY  GENERATION 

C         GIVE  INPUT  OF  FIRST  ELEMENT  NUMBER  OF  THE  LAYER, 
C         PROPERTY  ID  NUMBER,  THEN  INCREMENTAL  LAYER  NUMBER 
120  READ(1,2030)  NN,M,L 
N=NN 

IF(N.EQ.O)  GO  TO  140 
^.  -.     ■    ;         DO  9  J=1,NA 
I=J-1 

MPROP(N+NH*I)=M 
y^^--    "     •    9  CONTINUE 


Jit 
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':-^'  IF(L.EQ.O)     GO    TO    120 

IAUT0=1 
130   N=N+IAUT0 

DO    10    K=1,NA 
II=K-1 
10   MPROP(N+NH*II)=M 

IF(N.NE.L)  GO  TO  130 
GO  TO  120 
140  WRITE (2, 1040) 
1040  FORMAT(/, 3 X, 'MATERIAL  PROPERTIES  ARE  ASSIGNED',/) 
C  ****  STRESS-STRAIN  RELATION  INFORMATION 

150  READ(1,2040)  MP,C11(MP) ,C12(MP) ,C13(MP) ,C33(MP)  ,C44(MP) 
IF(MP.LT.MPTOTAL)  GO  TO  150 
DO  13  I=1,MPT0TAL 
MP=I 

CALL  DD(MP,C11,C12,C13,C33,C44,D) 
WRITE (2, 1071)  MP 

WRITE(2,1070)  ((D(K,J) ,J=1,4) ,K=1,4) 
13  CONTINUE 
1071  FORMAT (3X, 'STRESS-STRAIN  RELATIONS;  MATERIAL  ID=  ',15/) 
1070  FORMAT(4(F15.2,3X) ) 
2040  FORMAT(I5,5F12.4) 
C  ****  READ  DATA  OF  ULTIMATE  STRENGTH  OF  COMPOSITES        :- 
READ(1,2045)  XT, XC, XI , YT, YC, SI , S12 , S13 , S23 

2045  FORMAT (3F10. 3) 
WRITE(2,2046)  XT, XC,XI , YT, YC, SI , S12 , S13 , S23 

2046  FORMAT ( 3 X, 'ULTIMATE  STRENGTH  OF  THE  COMPOSITES ',/, 3X, 
r'     •'   -        +'XT=',F11.3,3X,  •XC=' ,F11.3,3X,  •XI=' ,F11.3,/,3X,  'YT=', 

+F11.3,3X, 'YC=',F11.3,3X, 'SI=',F11.3,/,3X, 'S12=' ,F10.3, 
+3X, 'S13=' ,F10.3,3X, 'S23=' ,F10.3,/) 
-      C  ****  CONTACT  INFORMATION 
READ(1,2050)  RB 
'■  2050  FORMAT  (FIO.  5) 

.  ■  >-  WRITE  (2, 1080)  RB 

1080  FORMAT ( 3 X, 'RADIUS' ,F8. 3, 'BALL  INDENTING  THE  PLATE',/) 
DO  I=1,NCAND 
NC0NT(I)=(NH+1)*(I-1)+1 
END  DO 

READ(1,2060)  NTRY 
-  2060  F0RMAT(I5) 

WRITE (2, 1091)  NTRY 
1091  FORMAT ( 3 X, 'MAXIMUM  CONTACT  PROGRESSION  NODE  ID. ',15,) 

READ(1,2060)  NREST  '- 

C  ****  INPUT  OF  DELAMINATION 

READ(1,2065)  NTCR  v.  [; 

.  '.'  IF(NTCR.EQ.O)  THEN  ;;"-' 

;^;;  WRITE  (2, 1044) 

1044  F0RMAT(3X, 'NO  DELAMINATION  PRESENT  FOR  THIS  RESTART',/) 
-'    •  GO  TO  172 

END  IF        . 
2065  FORMAT (15)  ' 

,     WRITE  (2, 1043)  NTCR  >        ;' 


T^^TT'-"  . 
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1043  FORMAT ( 3 X, 'NUMBER  OF  DELAMINATIONS ' , 15 ,/ , 3X, 

+ ' DELAMINATION  I . D. ', 3X, ' LENGTH (ELEMENT  N0.)',3X, 
+'STARTING  POINT(ELEMENT  I.D.)'//) 
DO  16  I=1,NTCR 

READ(1,2070)  NCRACK(I) ,  MECR(1,I) 
2070  F0RMAT(I5,I5) 

WRITE(2,1041)  I,  NCRACK(I) ,  MECR(1,I) 

1041  FORMAT(10X,I5,10X,I5,10X,I5) 

16  CONTINUE  '^i 

KSOLID=0  y 

DO  21  J=1,NTCR 

DO  12  I=2,NCRACK(J) 

MECR(I,J)=MECR(1,J)+NH*(I-1) 
12  CONTINUE 

DO  14  I=1,NCRACK(J) 

MNODE(MECR(I,J) ,2)=N0DES+I 

IF(I+1.LE.NCRACK(J) )  THEN 

MNODE(MECR(I,J) ,1)=N0DES+I+1 

END  IF 

RC(N0DES+I)=RC(MN0DE(MECR(I,J)+1,3) ) 

ZC (NODES+I) =ZC (MNODE (MECR (I , J) +1 , 3 ) ) 

NBC(N0DES+I,1)=1 

NBC (NODES+I, 2) =1 

NS0LID(I+KS0LID,1)=MN0DE(MECR(I,J) ,2) 

NSOLID(I+KSOLID,2)=MNODE(MECR(I,J)+l,3)  -;;,.' 

14  CONTINUE  "■■'{ 
NODES=NODES+NCRACK ( J) 

KSOLI D=KSOLI D+NCRACK ( J ) 

DO  15  I=1,NCRACK(J)-1 

MNODE (MECR ( I , J ) , 4 ) =MNODE ( MECR ( I , J ) +NH , 3 ) 

15  CONTINUE 

IF(RC(MN0DE(MECR(1,J)  ,2)  )  .EQ.0.0)  THEN  ,   •  -  i' 

NBC(MNODE(MECR(l,J) ,2) ,1)=0 

END  IF  _: 

21  CONTINUE  ■^^ 

DO  J=1,NTCR 

DO  I=1,NCRACK(J)  ', 

KDELA(MECR(I,J))=1  .      r 

END  DO 
END  DO 
»       WRITE(2,1042)  ( (NSOLID(I , J) , J=l , 2) , 1=1, KSOLID) 

1042  FORMAT (2 (15, 2X)) 
172  CONTINUE 

C  ****  INPUT  OF  FAILED  MATRIX  ELEMENTS 
READ(1,2066)  NRMF 
2066  F0RMAT(I5) 

IF(NRMF.EQ.O)  THEN 
WRITE (2, 1045) 
1045  F0RMAT(/,3X, 'NO  MATRIX  FAILURE  IN  THIS  RESTART',/) 

GO  TO  170  ■ 

END  IF  -^'^ 

WRITE (2, 1046)  NRMF  „ ,  ^ 
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1046  FORMAT ( 3 X, 'TOTAL  NO.  OF  STACKS  OF  FAILED  MATRIX', 15, 
+/f3X,'N0.  OF  ELEMENTS  IN  EACH  STACK' , 3X, ' FIRST  ELEMENT 
+ID.',/) 

DO  33  I=1,NRMF 
READ(1,2077)  NRMF1,NRMF2 
WRITE (2, 1047)  NRMF1,NRMF2 
2077  FORMAT(2I5) 

1047  FORMAT(25X,I5,10X,I5) 
DO  33  J=1,NRMF1 
DUM=NRMF2+ ( J-1) *NH 
MPR0P(DUM)=2 
KFAIL(DUM)=1 

33  CONTINUE 

170  CONTINUE 
C  ****  INPUT  OF  FIBER-FAILED  ELEMENTS 

READ (1,2066)  NRFF 
IF(NRFF.EQ.O)  THEN 
WRITE (2, 1051) 
1051  F0RMAT(/,3X, 'NO  FIBER  FAILURE  IN  THIS  RESTART',/) 
GO  TO  171 
END  IF 
WRITE (2, 1048)  NRFF 

1048  FORMAT (3X, 'TOTAL  NO.  OF  STACKS  OF  FAILED  FIBER' ,I5,/,3X, 
+'N0.  OF  ELEMENTS  IN  EACH  STACK' , 3X, ' FIRST  ELEMENT 
+ID.',/) 

DO  34  1=1, NRFF  "  , 

READ(1,2077)  NRFF1,NRFF2 
WRITE (2, 1049)  NRFF1,NRFF2 

1049  FORMAT (25X, 15, lOX,  15) 
DO  34  J=1,NRFF1 

DUM=NRFF2+(J-1)*NH  , 

MPR0P(DUM)=3  -    '# 

KFAIL(DUM)=3 

34  CONTINUE 

171  CONTINUE 

C  ****  NBC(N,I)  GIVES  GLOBAL  EQUATION  NUMBER  OF 

C      N;  NODE  NUMBER 

C      I;  D.O.F.  ID.  NUMBER  (R=l,  Z=2) 

DO  22  I=1,NCAND 

NBC(NCONT(I) ,2)=0 
22  CONTINUE 

NEQ=NCAND 

DO  2  5  N=l, NODES 

DO  25  1=1,2 

IF(NBC(N,I) .EQ.O)  GO  TO  25 

NEQ=NEQ+1 

NBC(N,I)=NEQ 
25  CONTINUE 

DO  28  I=1,NCAND  .  ,^ 

NBC  ( NCONT  ( I )  ,  2 )  =1  -Ji^, 

28  CONTINUE 
C  ****  PRINTING  GENERAL  INPUT  INFORMATION 


it-.c 
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WRITE (2, 1030) 
1030  FORMAT(//, 3 X, 'ELEMENT  DATA',//) 
DO  11  I=l,NEIiMS 

WRITE(2,1050)  I, (MNODE(I,J) ,J=1,4) ,MPROP(I) 
11  CONTINUE 
1050  FORMAT (3X, 'ELEMENT  ID' , 14 , 3X, 'NODES ' , 4 (3X, 13) , 3X, 
+'  MATERIAL  PROPERTY  ID ',14) 
DO  29  1=1, NODES 

WRITE(2,1100)  I,RC(I) ,ZC(I) ,NBC(I,1) ,NBC(I,2) 
29  CONTINUE 
1100  FORMAT  (3X,  'NODE  ID=',I4,3X,  'R-CORD' ,  3X,  F7  .  3  ,  3X,  'Z-CORD' 
+  ,3X,4F7.3,3X,  'R-GLOBAL  EQUATION  NUMBER  '  ,I4,3X,  'Z-G.E.N 
+',I4) 
C  ****  CALL  GLOBAL  STIFFNESS  CONSTRUCTION  AND  CONDENSATION 
CALL  GLOBS ( AX, BX,CX, CORD, RC,ZC,MPROP,NCAND, NBC, MNODE, 
+NEQ , NELMS , Cll , C12 , C13 , C3  3 , C44 , NSOLID , KSOLID , NTCR) 
C  ****  START  CONTACT  PROGRESSION 
KC0NT=1 

IF(NREST.NE.O)  THEN 
WRITE (2, 1201)  NREST 
1201  FORMAT ( 3 X, 'RESTART  PROGRAM  AT  KCONT=' , 15,/ , 3X, 
+' RESTARTING  CONTACTED  NODE  DISPLACEMENTS',/) 
KC0NT=NREST+1 
DO  K=1,KC0NT 
READ(1,2080)  DISP(K) 
I  '       ',!  '   2080  FORMAT  (F2  0.12)  i 

WRITE (2, 2080)  DISP(K) 
END  DO 
*■  ;■  END  IF 

WRITE  (2, 1900)  .^ 

-  999  CONTINUE  .  A 

IF(NREST.EQ.O)  THEN  ^ 

WRITE  (2, 1200)  KCONT,NCONT(KCONT)  "^^ 

WRITE (2, 1900)  H 

WRITE (2, 1900)  '  ■ 

END  IF 
888  NCR=NCAND-KCONT  ; 

IF (NREST. EQ.O)  THEN  ■   ' 

1200  FORMAT  (/,3X,  'PROGRESS  OF  THE  CONTACT' ,14,/,  3X,  'CONTACTED 

+AT  ',13, 'NODE  AND  NEXT  NODE  WILL  BE  CONTACTED')  ,   . 

END  IF  '         ' 

C  ****  CONDENSING  [AX]  INTO  [AA] , [AB] , [AC] 
CALL  CONCEN ( AX , AA , AB , AC , KCONT , NCAND) 
C  ****  [  AC  ] {  Zm  }=  -  [AB]^  {1,  l,...l}^ 

IF (NREST. NE.O)  THEN  <> 

DO  1=1, NCR 

DUM=0 . 0  ' 

DO  K=l, KCONT 

DUM=DUM-AB(K,I)*DISP(K)  j 

END  DO  ,i^ 

zM(i)=DUM  : 

END  DO  i 
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ELSE 

DO  30  1=1, NCR 

DUM=0.0 

DO  31  K=1,KC0NT 

DUM=DUM-AB(K,I) 
31  CONTINUE 

ZM(I)=DUM 
30  CONTINUE 

END  IF 

WRITE(2,1210)  (ZM(I) ,I=1,NCR) 
1210  FORMAT ( 3 X, 'INPUT  TO  SOLVE  ZM' ,/ , 5 (2X, F20 . 4) ) 
C  ****  SOLVE  LINEAR  EQUATION 
C       ZM  IS  TREATED  LIKE  INPUT  FORCE 
C       IMSL  PACKAGE  EMPLOYED 

IA=20 

IDGT=0 

MM=1 

CALL  LEQT1F(AC,MM,NCR,IA,ZM,IDGT,WKAR,IER) 

IF  (NREST.EQ.O)  THEN 

WRITE(2,1220)  (ZM(I) , I=1,NCR) 

1220  FORMAT (3X, 'ZM:  DISP.  TO  THE  UNIT  INPUT ' ,/, 5 (2X, F15 . 8) ) 
ELSE 

WRITE(2,1221)  (ZM(I) ,I=1,NCR) 

1221  FORMAT (3X, 'RESTARTED  DISPLACEMENTS  OF  NODES  NOT 
+CONTACTED  YET ' ,/ , 3X, 6 (F12 . 6, 2X) ) 

'  DELTA=ZM(1)-DISP(KC0NT) 

GAP=DSQRT(RB**2-RC(NC0NT(KC0NT) ) **2) -DSQRT (RB**2- 
+RC ( NCONT ( KCONT+ 1 ) ) *  *  2 ) 

IF( (DELTA+GAP) .LE.0.0)  THEN 
KC0NT=KC0NT+1 
DISP(KC0NT)=ZM(1) 

WRITE(2,1211)  DELTA,  GAP, DISP (KCONT) , NCONT (KCONT) 
(.  1211  FORMAT  ( 2  X,  'AUTOMATIC  NEW  CONTACT  BETWEEN  INDENTER 

^>  +AND  PLATE' ,/,3X,  'Dj-Di' ,F10.5,  'GAP' ,F10.5,/,3X, 

+ 'DISP (KCONT) =' ,F12. 6, 3X,  'NEW  CONTACTED  NODE ', 15) 
GO  TO  888 
END  IF 
END  IF 
'- ,       C  NOW,  Zm  IS  OUTPUT  DISPLACEMENTS 

C  [AC]  IS  (NCR  X  NCR)  DIMENSION 

C  (Zm)   IS  (NCR  X  1)  DIMENSION 

NR=NEQ-NCAND 
,  -         IF(NREST.NE.O)  THEN 
DO  1=1, KCONT 
'■  DUM1=0.0 

DUM2=0.0 
.A-*>  ':  '   ^  ,  DO  J=l, KCONT 

'^'    '"'■■■■       •^':'        DUM1=DUM1+AA(I,J)*DISP(J) 
END  DO 
>  ,  \  •     '  DO  K=1,NCR 

DUM2=DUM2+AB ( I , K) *ZM (K) 
END  DO 

""  ,-  .--  ■'  >  .«    -         '  .i  ■  ' 

.■■  '  ■ "'  ■     ^v    ."  „..i   ■': 


*•-'». 
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-  F(I)=DUM1+DUM2 

END  DO  ■■/■:■ 

"•.  DO  JJ=1,NCR 

■;  ■:;  DISP(KCONT+JJ)=ZM(JJ) 

"'  END  DO  ' 

DUM=0 .0 
DO  J=l , KCONT 
DUM=DUM+F(J) 
END  DO 

WRITE (2, 1900) 
WRITE (2, 1231)  DUM 
WRITE(2,1232)  (F (I) , 1=1, KCONT) 
1231  FORMAT (3X, 'RECALCULATED  TOTAL  FORCE  AFTER  FAILURE', 
■^-  +F10.3,/,  3 X,  'NODAL  FORCES' ,/) 

L^         1232  FORMAT(3X,F20.8)  :\ 

DO  1=1, NR 
;,  .  DUM=0.0 

!  DO  K=1,NCAND 

DUM=DUM+BX (K, I) *DISP (K) 
END  DO 
;  DISP(I+NCAND)=DUM 

END  DO 
:^       C  ****  PRINTING  RECALCULATED  DISPLACEMENTS  ' 

DO  IJ=2,NH+1 

DUM=DISP(NBC(IJ,2) )-DISP(l) 
WRITE(2,1234)  ZC (IJ) , DISP(NBC(IJ, 2) ) , DUM 
END  DO  ^ 

WRITE (2, 1235)  RC (NCONT (KCONT+1) ) 
WRITE (2, 1900) 
,  1234  F0RMAT(3X, 'Z-CORD' ,F8. 4, 2X, 'DISP' ,F12. 8, 2X, 'RELATIVE 
+INDENTATION ' , F12 . 8 ) 
1235  FORMAT (/,3X, 'RECALCULATED  RADIUS' ,F8. 5,/) 
C  ****  CHECK  DAMAGES (CALCULATE  ELEMENT  STRESSES  AND  APPLY 
C        MODIFIED  HASHIN'S  CRITERIA) 
-;  .  CALL  RECOVER(CORD,RC,ZC,MNODE,MPROP,Cll,C12,C13,C33,C44, 

' "  '  +D, NBC, DISP, STRESS, STRAIN, NA,NH, GAUSS, XT, YT,XC,YC,S 12, 

,;f.  +S13,S23,XI,SI,SA,NELMS,KDELA,KFAIL) 

NREST=0 
GO  TO  999 
END  IF 
C  ****  CLOSE  GAP  AT  THE  NEXT  CONTACT  CANDIDATE  NODE 

GAP=DSQRT(RB**2-RC (NCONT (KCONT) ) **2) -DSQRT(RB**2- 
•"  - '  +RC  (  NCONT  (  KCONT+ 1 )  )  *  *  2 ) 

CM= (DISP (KCONT+1) -DISP (KCONT) +GAP) /( 1. -ZM ( 1) ) 
WRITE (2, 1230)  CM 
1230  FORMAT(/,3X, '  ADJUSTMENT  FACTOR  TO  CLOSE  GAP',F20.8) 
IF(CM.LE.O.O)  THEN 
WRITE (2, 1233) 
■  1233  FORMAT ( 3 X, 'WARNING;  CONTACT  PROCEEDS  WITHOUT 
+EXERTING  ANY  FORCE') 
v        END  IF 
l.„   ■•  ."    c  ****  CALCULATE  DISPLACEMENTS  AND  FORCES 
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DO  35  I=1,KC0NT 
ZU(I)=CM 

35  CONTINUE 

DO  36  J=1,NCR 

ZU (KCONT+J) =ZM (J) *CM 

36  CONTINUE 

WRITE(2,1240)  (ZU(I) , I=1,NCAND) 
1240  F0RMAT(/,3X, 'UPPER  SURFACE  INCREMENTAL  DISPLACEMENT',/, 
+5(F12.8,2X) ) 
C  ****  {RF}  =  CM*([AA]{1,1,  .  .  .1)^  +  [AB]{Zin)) 
DO  40  I=1,KC0NT 
-.  DUM=0.0 

DO  41  K=1,KC0NT 
DUM=DUM+AA(I,K) 

41  CONTINUE 
DUM1=0 . 0 

DO  42  J=1,NCR 
DUM1=AB(I,J)*ZM(J)+DUM1 

42  CONTINUE 

RF ( I ) = ( DUM+DUMl ) *CM 
40  CONTINUE 

WRITE(2,1250)  (RF (I) , 1=1 ,NCAND) 
1250  FORMAT (/, 3 X, 'INCREMENTAL  NODAL  FORCE ' ,/ , 5 (IX, F20 . 4) ) 
C  ****   (  Zr  }   =   [  X  ]  {  Zu  },  SEE  SUBROUTINE  GLOB 
DO  45  1=1, NR 
DUM=0 . 0 

DO  46  K=1,NCAND 
DUM=DtJM+BX  (K,  I )  *ZU  (K) 
46  CONTINUE 

ZR(I)=DUM 
45  CONTINUE 
C  ****   TOTAL  DISPLACEMENTS  AND  NODAL  FORCES 
DO  55  I=1,NCAND 
DISP(I)=DISP(I)+ZU(I) 
F(I)=F(I)+RF(I) 
55  CONTINUE 

DO  60  1=1, NR 
J=I+NCAND 

DISP(J)=DISP(J)+ZR(I) 
60  CONTINUE 

WRITE(2,1280)  (F(I) , I=1,NCAND) 
1280  FORMAT (/,3X,  ' UPDATED  CONTACT  NODAL  FORCE  '  ,/ ,  5  (IX,  F15.  4)  ) 

WRITE(2,1290)   (DISP (I) , I=1,NCAND) 
1290  FORMAT(/, 3 X, 'UPDATED  DISPLACEMENTS ',/, 5 (IX, F15. 8) ) 
C  ****  CALCULATE  ELEMENT  STRAINS  AND  STRESSES 

CALL  RECOVER (CORD, RC,ZC,MNODE,MPROP,Cll,C12,C13,C33, 
+C44,D, NBC, DISP, STRESS, STRAIN, NA,NH, GAUSS, XT, YT,XC,YC, 
+S12 , S13 , S23 , XI , SI , SA, NELMS , KDELA, KFAIL) 
C  ****  NOTATION  OF  STRESS  DATA,  STRESS (N,L) 
C         N;  NODE  IDENTIFICATION 
C  STRESS (N,l)  ;  o^ 

C  STRESS (N, 2)  ;  o. 
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C  STRESS  (N,  3)  ;  x ^^ 

C  STRESS (N, 4)  ;  Og 

C  ****  GENERAL  CONTACT  INFORMATION 
DXJM=0 .  0 

DO  80  I=1,KC0NT 
DUM=F(I)+DXJM 
80  CONTINUE 

FXX(KCONT)=DUM 
WRITE (2, 1902) 
1902  FORMAT ( 3X, *****  CONTACT  INFORMATION  •) 

WRITE(2,1400)  NC0NT(KC0NT+1) ,RC(NC0NT(KC0NT+1) ) ,DISP(1) , 
+FXX(KCONT) 
1400  F0RMAT(4X,  'N0DE',5X,  'CONTACT  RADIUS', lOX,  'TOP  DEFLECT', 
+3X, 'CONTACT  FORCE' ,/ , 4X, 14 , 2X, F20 . 5, IX, F20 . 8 , 3X, F30 . 4) 
DUM=DISP(1)-DISP(NBC(1+NH,2) ) 
WRITE (2, 1901)  DISP(NBC(1+NH,2) ) , DUM 
1901  F0RMAT(3X,  '  BOTTOM  DEFLECTION= ', F2  0 . 12 , 3X,  ' INDENTATION= ' 
+,F20.12) 
WRITE (2, 1900) 
DO  IJ=2,NH+1 

DUM=DISP(NBC(IJ,2) )-DISP(l) 
WRITE(2,1234)  ZC(IJ) , DISP(NBC (IJ, 2) ) , DUM 
END  DO 

WRITE (2, 1900) 
1900  FORMAT(3X, '*************************************•) 

WRITE (2, 1500) 
1500  FORMAT (/,3X,  '  RADIUS  ',5X,  'AREA' ,10X,  'Szz' ,10X,  'THEORY') 
WRITE (2, 1900) 
PI=3.141592625D0 
DO  85  I=1,KC0NT 
J=(I-1)*NH+1 
IF(I.EQ.l)  THEN 

AREA=RC(NC0NT(2) )**2*PI/4. 
ELSE 

AREA=RC(NCONT(I) )*PI*2.*XAC 
END  IF 

PRESS=STRESS ( J , 2 ) 
DUM=RC(NC0NT(KC0NT+1) ) **2 

THE0=-1 . 5*FXX (KCONT) /PI/DUM*SQRT ( 1 . -RC (NCONT (I) ) **2/DUM) 
WRITE(2,1550)  RC(NCONT(I) ) , AREA, PRESS, THEO 
1550  FORMAT(3X,F8.4,F9.4,F20.4,F20.4) 
85  CONTINUE 
C  ****  ITERATE  TO  NEXT  CONTACT  NODE  UNTIL  INDENTER  CONTACTS 
C       NCONT (NTRY)  NODE 
WRITE (2, 1900) 
KC0NT=KC0NT+1 

IF (KCONT. LE. NTRY)  GO  TO  999 
STOP 
END 
C 

SUBROUTINE  AXISYM(CORD, S , D) 
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IMPLICIT  REAL*8 ( A-H , 0-Z ) 

DIMENSION  S(8,8) ,CORD(2,4) ,D(4,4) ,B(4,8) ,DB(4) 

DIMENSION  GAUSS (4) , WEIGHT (4) 
C  ****  GAUSS  -LEGENDRE  INTEGRATION  COEFFICIENTS 

DATA  WEIGHT/  1. ODO , 1 . ODO , 0 . , 0 ./ 

DATA  GAUSS/  0. 5773502692D0, -0. 5773502692D0 , 0 . , 0 ./ 
C  ****  INITIATE  STIFFNESS  MATRIX 

DO  10  1=1,8 

DO  10  J=l,8 

S(I,J)=0.0 
10  CONTINUE 
C  ****  EVALUATE   MATRIX  B 

PI=3.141592625D0 

DO  20  1=1,2 

R=GAUSS(I) 

DO  20  J=l,2 

Z=GAUSS(J) 

CALL  BJ (CORD, B,DET,R,Z, RADIUS, D) 

WT=WEIGHT(I) *WEIGHT(J) *RADIUS*DET*2 . 0*PI 

DO  30  JJ=1,8 

DO  40  K=l,4 

DB(K)=0.0 

DO  40  L=l,4 
40  DB(K)=DB(K)+D(K,L)*B(L,JJ) 

DO  50  M=JJ,8 

STIFF=0.0 

DO  60  LL=1,4 
60  STIFF=STIFF+B(LL,M)*DB(LL) 
50  S(M,JJ)=S(M,JJ)+STIFF*WT 
30  CONTINUE 
2  0  CONTINUE 

DO  90  J=l,8 

DO  90  I=J,8 

S(J,I)=S(I,J) 
90  CONTINUE 

RETURN 

END 
C 

SUBROUTINE  BJ (CORD, B, DET,R, Z , RADIUS , D) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSIONCORD(2,4)  ,B(4,8)  ,H(4) ,P(2,4) ,XJI(2,2) ,XJ(2,2) 

RP=1.0+R 

SP=1.0+Z 

RM=1.0-R 

SM=1.0-Z 

C  ****  INTERPOLATION  FUNCTION 
H(1)=0.25*RP*SP 
H(2)=0.25*RM*SP 
H(3)=0.25*RM*SM 
H(4)=0.25*RP*SM 
P(1,1)=0.25*SP 
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P(l,2)=-P(l,l) 

P(1,3)=-0.25*SM 

P(l,4)=-P(l,3) 

P(2,1)=0.25*RP 

P(2,2)=0.25*RM 

P(2,3)=-P(2,2) 

P(2,4)=-P(2,l) 

DO  30  1=1,2 

DO  30  J=l,2 

DUM=0 . 0 

DO  20  K=l,4 
20  DUM=DUM+P(I,K)*CORD(J,K) 
30  XJ(I,J)=DUM 

DET=XJ(1,1)*XJ(2,2)-XJ(2,1)*XJ(1,2) 

IF(DET.GT. 0.0000001)  GO  TO  40 

WRITE (2, 1000)  ((CORD (I, J) ,J=1,4) ,1=1,2) 
1000  F0RMAT(3X, 'JACOBIAN  FAILED  AT  ' ,/ , 3X, 8 (F9 . 4 , IX) ) 

STOP 
40  DUM=1./DET 

XJI (1,1) =XJ (2,2) *DUM 

XJI(1,2)=-XJ(1,2)*DUM 

XJI  (2  , 1)  =-XJ  (2  , 1)  *DU]yi 

XJI(2,2)=XJ(1,1)*  DUM 
C  ****  EVALUATE  GLOBAL  DERIVATIVE  OPERATOR 

KK=0 

DO  60  K=l,4 

KK=KK+2 

B(1,KK-1)=0. 

B(1,KK)=0. 

B(2,KK-1)=0. 

B(2,KK)=0. 

DO  50  1=1,2 

B(1,KK-1)=B(1,KK-1)+XJI(1,I)*P(I,K) 
50  B(2,KK)=B(2,KK)+XJI(2,I)*P(I,K) 

B(3,KK)=B(1,KK-1) 
60  B(3,KK-1)=B(2,KK) 

RADIUS=0.0 

DO  70  K=l,4 
70  RADIUS=RADIUS+H(K)*C0RD(1,K) 

DUM=1. /RADIUS 

KK=0 

DO  90  K=l,4 

KK=KK+2 

B(4,KK)=0. 
90  B(4,KK-1)=H(K)*DUM 

RETURN 

END 
C 

SUBROUTINE  DD (MP, Cll , C12 , C13 , C3  3 , C44 , D) 

IMPLICIT  REAL*8(A-H,0-Z) 
C  ****  STRESS-STRAIN  RELATIONS 

DIMENSION  D(4,4) 


^ 


■  127 

:-'-■■■,:,  I'- 

V   "       DIMENSION  Cll(l) ,C12(1) ,C13(1) ,C33(1) ,C44(1) 
^-      ■    DO  10  1=1,4  ,W.. 

DO  10  J=l,4     '■    ,  ', 

D(I,J)=0.0 
■■*!..      -;':    10  CONTINUE  "^  •  - 

>:.^  ■  D(1,1)=C11(MP) 

D(1,2)=C13(MP)  ^   1 

D(1,4)=C12(MP) 
"    '        D(2,2)=C33(MP) 
D(2,4)=C13(MP) 
D(3,3)=C44(MP) 
'   V   .'  '     D(4,4)=C11(MP) 
'  • ;        DO  20  1=1,4 
DO  20  J=I,4 
D(J,I)=D(I,J) 
20  CONTINUE 

•  •'        RETURN 

END 

SUBROUTINE  CLEAR (ZZ,N,M )  .  ! 

.■     IMPLICIT  REAL*8(A-H,0-Z)  ! 

DIMENSION  ZZ(1,1)  ,'    j. 

f       DO  1  1=1,  N  ,y'.  '"  ■;.  . 

^  '   DO  1  J=1,M  ':'%  '  '  «-  :         ■:■  .  "    V  " 

J       ZZ(I,J)=0.0  ,  ^: 

^   ^' -     1  CONTINUE       'i        ■  ^..    .:.-.  : 

'.■:\     .  '        RETURN    i(,  '  '     "  •'  ■'■ 

*  END       '\.'  •:.      ■ 
C 

SUBROUTINE  GLOBS (AX, BX, CX, CORD, RC, ZC,MPROP,NCAND, NBC, 
+MNODE ,  NEQ ,  NELMS ,  Cll ,  C12  ,  C13  ,  C33 ,  C44  , NSOLID,  KSOLID,  NTCR) 
IMPLICIT  REAL*8(A-H,0-Z) 
-DIMENSION  AX(20,20) ,BX(20,1250) ,CX(1250,1250) ,D(4,4) 
'  DIMENSION  C0RD(2, 4) ,RC(650) ,ZC(650) ,MPROP(550) 

DIMENSIONNBC(650,2) ,MNODE (550, 4) ,X(1250,20) ,WKAR(5500) 
DIMENSION  Cll(l) ,C12(1) ,C13(1) ,C33(1) ,C44(1) ,S(8,8) 
DIMENSION  NSOLID (150, 2) 
C  ****  CONSTRUCTING   GLOBAL  STIFFNESS 
A'  .  NR=NEQ-NCAND 

/-:.  WRITE  (2, 1000)  NR 

.  •'  •.   ,  1000  FORMAT(/,  3  X,  'SUBROUTINE  GLOBS  IS  CALLED  IN  ,NR=',I5) 
CALL  CLEAR (AX, NCAND,NCAND) 
CALL  CLEAR(BX,NCAND,NR) 
;.*:'  -^   ■;.-•    CALL  CLEAR (CX,NR,NR) 
"^        IC=0 

DO  10  N=l, NELMS 
DO  20  JJ=1,4 

C0RD(1, JJ)=RC(MNODE(N, JJ) ) 
\        C0RD(2,JJ)=ZC(MN0DE(N,JJ) ) 
20  CONTINUE 
C  ****  CALL  SUBROUTINE  DD  WHEN  MATERIAL  PROPERTY  OF  A  ELEMENT 
C,^.  ,    .   HAS  CHANGED 


V*.;. 
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MP=MPROP(N) 

IF(IC.EQ.MP)  GO  TO  100 
CALL  DD(MP,C11,C12,C13,C33,C44,D) 
100  IC=MP 

CALL  AXISYM(CORD,S,D) 
DO  10  1=1,8 
K=(I+l)/2 
L=I-2*K+2 

NROW=NBC (MNODE (N, K) , L) 
IF(NROW.EQ.O)  GO  TO  10 
DO  10  J=l,8 
K=(J+l)/2 
L=J-2*K+2 

NCOL=NBC (MNODE (N,K) , L) 
IF(NCOL.EQ.O)  GO  TO  10 
IR=NROW-NCAND 

IC=NCOL-NCAND  .{ 

IF(IR.GT.O.AND.IC.GT.O)THEN 
*        CX(IR,IC)=CX(IR,IC)+S(I,J) 

ELSE  IF(IR.LE.O.AND.IC.GT.O)THEN 
BX (NROW, IC) =BX (NROW, IC) +S ( I , J) 
ELSE  IF(IR.LE.O.AND.IC.LE.O)THEN 
AX (NROW, NCOL) =AX (NROW, NCOL) +S ( I , J) 
■:  ELSE  V   f 

GO  TO  10 
END  IF 
10  CONTINUE 
C  ****  IMPLEMENT  DELAMINATIONS 

IF(NTCR.EQ.O)  THEN  V^, 

GO  TO  200 
END  IF 

XK=9000000000.0 
DO  15  I=1,KS0LID 
IR=NBC (NSOLID ( I , 1 ) , 2 ) -NCAND 
IC=NBC (NSOLID ( I , 2 ) , 2 ) -NCAND 
CX ( IR, IR) =CX ( IR , IR) +XK 
CX(IC,IC)=CX(IC,IC)+XK 
CX ( IR, IC) =CX (IR, IC) -XK 
CX ( IC , IR) =CX ( IC , IR) -XK 
15  CONTINUE 
200  CONTINUE 
C  ****  CONDENSE;  [K]*{Zu)  ={RF} 
C         [X]  =  -  [CX]'^  [BX]^ 
C         [X]  HAS  DIMENSION  OF  (NR  X  NCAND) 
C  {Zr}=[X]{Zu} 

C  [AA] {Zu}={RF},  WHERE 

C  [AA]  =  [AX]-[BX]  [CX]"''[BX]^ 

C  ****  FIRST,  DEFINE  [X]  AS  -[BX]^ 
DO  30  1=1, NR 
DO  30  J=l, NCAND 
X(I,J)=-BX(J,I) 
30  CONTINUE 
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C  ****  NOW,  MAKE  [X]=-[CX]'''[BX]^ 

IDGT=0 

IA=1250 

CALL  LEQT1F(CX,NCAND,NR,IA,X,IDGT,WKAR,IER) 
C  ****  MAKE  AA  MATRIX  (CONDENSED  FORM  ) 

DO  40  I=1,NCAND 

DO  40  J=1,NCAND 

DUM=0 . 0 

DO  50  K=1,NR 

DUM=DUM-BX (I , K) *X (K, J) 
50  CONTINUE 

AX ( I , J ) =AX ( I , J ) - DUM 
40  CONTINUE 

DO  60  I=1,NCAND 

DO  60  J=1,NR 

BX(I,J)=X(J,I) 
60  CONTINUE 
C  ****  END  OF  CONDENSATION 

RETURN 

END 

SUBROUTINE  CONCEN(AX, AA, AB, AC,KCONT,NCAND) 

IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  AX(20, 20), AA(20, 20) 

DIMENSION  AB(20,20) ,AC(20,20) 

NCR=NCAND-KCONT 

CALL  CLEAR  (AA,  KCONT,  KCONT)        / /' - 

CALL  CLEAR (AB,KCONT, NCR) 

CALL  CLEAR (AC, NCR, NCR)  ,  - 

DO  10  I=1,NCAND 

DO  10  J=1,NCAND 

IF ( I . LE . KCONT . AND . J . LE . KCONT) THEN 

AA(I,J)=AX(I,J) 
ELSE  IF ( I. LE. KCONT. AND. J. GT. KCONT) THEN 

AB ( I , J-KCONT) =AX ( I , J) 
ELSE  IF ( I. GT. KCONT. AND. J. GT. KCONT)  THEN 

AC (I -KCONT, J-KCONT) =AX (I , J) 
ELSE 

GO  TO  10 
END  IF 
10  CONTINUE 
RETURN 
END 
C 

SUBROUTINE  RECOVER ( CORD, RC, ZC,MNODE,MPROP, CI 1 , C12 , 
+C13 , C33 , C44 , D, NBC, DISP, STRESS , STRAIN, NA, NH, GAUSS , 
+XT,YT,XC,YC,S12,S13,S2  3,XI,SI,SA,NELMS,KDELA,KFAIL) 
IMPLICIT  REAL*8(A-H,0-Z) 

DIMENSION  CORD(2,4) ,RC(650) ,ZC(650) ,MNODE (550 , 4) , 
DIMENSION  Cll(5) ,C12(5) ,C13(5) ,C33(5) ,C44(5) ,D(4,4) 
DIMENSIONDISP (12 50 ), STRESS (550,4) ,STRAIN( 550,4) , QUE (8) 
DIMENSION  ST(4) ,SA(4) ,SN(4) ,ASN(4) ,GAUSS(4) ,B(4,8) 
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DIMENSION  NBC(650,2) ,KFAIL(550) ,KDELA(550) ,MPROP(550) 

IC=0 

DO  10  K=1,NA 

DO  20  1=1, NH 

II=(K-1)*NH+I 

DO  30  J=l,4 

C0RD(1,J)=RC(MN0DE(II,J)) 

C0RD(2,J)=ZC(MN0DE(II,J)) 
30  CONTINUE 

MP=MPROP(II) 

IF(IC.EQ.MP)  GO  TO  1000 

CALL  DD (MP, Cll , C12 , C13 , C33 , C44 , D) 
1000  IC=MP 
C  ****  CALCULATE  ELEMENT  NODAL  DISPLACEMENT 

DO  40  L=l,8 

MM=(L+l)/2 

NN=L-2*MM+2 

IF(NBC(MNODE(II,MM) ,NN) .EQ.O)  THEN 

QUE (L) =0.0 

ELSE 

QUE(L)=DISP(NBC(MNODE(II,MM) ,NN) ) 

END  IF 
40  CONTINUE 

DO  45  J=l,4 

SA(J)=0.0 
45  CONTINUE  ■  .  ■;  •  ,^. 

DO  50  M=l,2  "^  ' 

DO  50  N=l,2 

R=GAUSS(M) 

Z=GAUSS(N) 

CALL  BJ( CORD, B,DET,R,Z, RADIUS, D) 
C  ****  CALCULATE  ELEMENT  STRAIN 

DO  60  IN=1,4 

DUM=0.0 

DO  65  IM=1,8 

DUM=DUM+B(IN, IM) *QUE (IM) 
65  CONTINUE 

SN(IN)=DUM 
60  CONTINUE 
C  ****  CALCULATE  STRESS 

DO  70  LL=1,4 

DUM=0 . 0 

DO  75  JJ=1,4 

DUM=DUM+D(LL, JJ) *SN(JJ) 
75  CONTINUE  ; 

ST(LL)=DUM 
70  CONTINUE 

DO  80  J=l,4 

SA(J)=SA(J)+ST(J) 

ASN ( J) =ASN ( J) +SN ( J) 
80  CONTINUE 
50  CONTINUE  -', 
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DO  85  J=l,4 
SA(J)=SA(J)*0.25 
ASN(J)=ASN(J)*0.25 
STRESS(II,J)=SA(J) 
STRAIN (II , J) =ASN ( J) 
85  CONTINUE 

CALL  CRIT(II,XT,YT,XC,YC,S12,S13,S23,XI,SI,SA, 
+KFAIL,KDELA,NELMS) 
20  CONTINUE 
10  CONTINUE 
DO  K=1,NA 

WRITE(2,1200)  (K-1)*NH+1, (KFAIL(I+ (K-1) *NH) ,I=1,NH) , 
+(KDELA(I+(K-1)*NH) ,I=1,NH) 
END  DO 
1200  FORMAT (IX, 13, 3X, 16 (IX, II) , 4X, 16 (IX, II) ) 
RETURN 
END 
C 

SUBROUTINE  CRIT (II , XT, YT,XC, YC, S12 , S13 , S23 ,XI , SI , SA, 
+KFAIL , KDELA , NELMS ) 
IMPLICIT  REAL*8 ( A-H , 0-Z ) 

DIMENSION  SA(4) ,IFFAIL(72) ,IMFAIL(72) ,KFAIL(550) , 
+KDELA(550) 
PI=3.14159262  5D0 
DO  5  1=1,72 

IFFAIL(I)=0  .,.  ;  , 

IMFAIL(I)=0 
5  CONTINUE 

DO  10  1=1,72 
XI=I 

X=PI/36.0*XI 
C 

C=COS(X)**2*SA(l)+SIN(X)**2*SA(4) 
D=(SIN(X) *COS (X) * (SA(4) -SA(1) )/S12) **2 
E=(COS(X)*SA(3)/S13)**2 

F=C+SA(2)  'i-    : 

G=(SIN(X)*SA(3) )**2-SIN(X)**2*SA(l)*SA(2)- 
+COS(X)**2*SA(4)*SA(2) 
C  ****  CHECKING  FIBER  FAILURE 
IF(C.LT.O.O)THEN 

IF(((C/XC)**2-1.0) .GE.0.0)  THEN 
IFFAIL(I)=1 
END  IF 
ELSE 

IF( ( (C/XT) **2+D+E) .GE.1.0)THEN 

IFFAIL(I)=1 
END  IF 
END  IF 
C  ****  CHECKING  MATRIX  FAILURE 
IF(F.GE.O.O)THEN 
H=(F/YT) **2+G/(S23**2)+D+E 
IF(H.GE.1.0)THEN 
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IMFAIL(I)=1 

...  END  IF 

■W  '  .       ELSE 

A=F*(yC/(2.*S23)**2-l./YC)  +  (F/(2.*S23))**2+G/S23**2+EH-E 

IF (A. GE. 1.0) THEN 
'-    /      IMFAIL(I)=1 
END  IF 
END  IF 
(,:'  10  CONTINUE 

C  ****  MATRIX  FAILURE  PRINTOUT.  PARTIAL  FAILURES  ARE  PRINTED 
C       HERE.  FOR  OTHER  INFORMATION,  PRINT  KFAIL(I) 
DO  50  1=1,72 

IF(IMFAIL(I) .EQ.l)  THEN 

MI=I 
GO  TO  110 
END  IF 
50  CONTINUE 
GO  TO  998 
110  CONTINUE 

IF (MI. EQ.l)  THEN 
DO  55  J=2,72 
IF(IMFAIL(J) .EQ.O)  THEN 
WRITE(2,1500)  II,J*5 
1500  FORMAT ( 3X, '**  MATRIX  FAILURE  AT ', 15, 3X, 'ELEMENT' , 3X, 
+  'ZERO  THROUGH' ,15,3 X,  'DEGREE  ,  SYMMETRIC) 
KFAIL(II)=2 
GO  TO  998 
END  IF 
55  CONTINUE 

KFAIL(II)=1  ,   . 

ELSE 

DO  60  J=MI+1,72 
IF(IMFAIL(J) .EQ.O)  THEN 
WRITE(2,1700)  II  ,MI*5,J*5 
1700  F0RMAT(3X, '**  MATRIX  FAILURE  AT  ', 15, 3X, 'ELEMENT ', 3X, 
+' FROM ',15,'  DEGREE  TO ',15,'  DEGREE  , SYMMETRIC  ') 
KFAIL(II)=2 
GO  TO  998 
END  IF 
60  CONTINUE 
END  IF 
998  CONTINUE 
C  ****  FIBER  FAILURE  PRINTOUT.  PARTIAL  FAILURES  ARE  PRINTED 
C       HERE.  FOR  OTHER  INFORMATION,  PRINT  KFAIL(I) 
DO  30  1=1,72 

IF(IFFAIL(I) .EQ.l)  THEN 

MI=I 
GO  TO  100 
END  IF 
30  CONTINUE 
;:  GO  TO  999 

,■'  100  CONTINUE  .  -; 
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IF(MI.EQ.l)  THEN 
DO  35  J=2,72 
IF(IFFAIL(J) .EQ.O)  THEN 
WRITE(2,1000)  II,J*5 
1000  FORMAT ( 3X, '$$$  FIBER  FAILURE  AT ', 15 , 3X, ' ELEMENT ', 3X, 
+ • ZERO  THROUGH • , 15 , 3X, ' DEGREE  ,  SYMMETRIC ' ) 
KFAIL(II)=4 
GO  TO  999 
END  IF 
35  CONTINUE 

KFAIL(II)=3 
ELSE 

DO  40  J=MI+1,72 
IF(IFFAIL(J) -EQ.O)  THEN 
WRITE(2,1200)  II  ,MI*5,J*5 
1200  FORMAT ( 3X, '$$$  FIBER  FAILURE  AT ', 15 , 3X, ' ELEMENT ', 3X, 
+' FROM ',15,'  DEGREE  TO ',15,'  DEGREE  , SYMMETRIC  •) 
KFAIL(II)=4 

GO  TO  999  V 

END  IF  ' ^ 

40  CONTINUE 
END  IF 
999  CONTINUE 
C  ****   DELAMINATION  PRINTOUT.  FOR  THE  INFORMATION,  PRINT 
C         KDELA(I),  I  =  ELEMENT  IDENTIFICATION 
IF(SA(2) .GT.0.0)  THEN 
B=(SA(2)/XI)**2+(SA(3)/SI)**2 
IF(B.GE.l.O)  THEN 
KDELA(II)=2 
END  IF 
ELSE 
B=(SA(3)/SI)**2 

IF(B.GE.l.O)  THEN 
KDELA(II)=1 
END  IF 
END  IF 
RETURN 
END 
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